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ABSTRACT 


Nonlinear effects are introduced in the dynamics of large space truss 
structures by the connecting joints which are designed with rather 
important tolerances to facilitate the assembly of the structures in 
space. The purpose of this work was to develop means to investigate the 
nonlinear dynamics of the structures, and particularly the limit cycles 
that might occur when active control is applied to the structures. An 
analytical method was sought and derived to predict the occurrence of 
limit cycles euid to determine their stability. This method is mainly 
based on the quasi- linearization of every joint using describing 
functions. This approach was proven successful when simple dynamical 
systems were tested. Its applicability to larger systems depends however 
on the amount of computations it requires, and estimates of the 
computational task tend to indicate that the number of individual 
sources of nonlinearity should be limited. Alternate einalytical 
approaches, which do not account for every single nonlinearity, or the 
simulation of a simplified model of the dynamical system should 
therefore be investigated to determine a more effective way to predict 
limit cycles in large dynamical systems with an important number of 
distributed nonlinearities. 
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CHAPTER ONE 
INTRODOCTION 

1.1 Presentation of some New Classes of Problems Associated with the 
Design of Large Space Structures. 

A new stage in space development is approaching that sees an 
important increase in space-based activities. This increase should 
result in spacecraft of much bigger dimension, where size will be mainly 
determined by the need to collect and transmit more radiative energy 
like solar power or radio signals, the only form of energy that can be 
exchanged in outer space. The Hoop/Column antenna, as well as the 
current Large Space Station configuration should have, for example, 
dimensions exceeding 100 meters, and comi>arable, or even bigger 
spacecraft, stations or satellites are expected in a more distaint 
future. 

Truss structures offer a satisfactory answer to the problem of 
sp>anning large distances or large areas, and they are extensively used 
in the design of new spacecraft: they constitute, for example, the 

backbone of the space station, or so-called ’'power tower" configuration. 
The reason for this is that they have a very high stiffness to mass 
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ratio, which is even more important in the construction of edifices in 
space than it is on Earth. 

However, the use of truss structures in space presents specific 
problems of vibration atnd instability. These problems arise from the 
inherent flexibility of structures which have very large geometrical 
characteristics, and which are built with a minimal amount of material, 
as weight is a limiting factor in the launch process. The situation is 
aggravated by the way the structures are assembled. Since it is not 
feasible to weld or bolt the different parts together in space, or to 
launch a structure in one piece, joints with rather large tolerances 
must be used so that structures can be erected by astronauts, or 
automatically deployed. The resulting stiffness is therefore reduced, 
and low frequency vibrations occur that may be very detrimental both for 
the hardware and for the completion of the mission. The behavior of the 
joints can even become predominant in the dynamics, and a lot of the 
benefits expected from a truss can be lost in joint dominated 
structures. A sometimes strongly nonlinear behavior of the joints 
induces other unwanted effects and complicates further the study of 
these problems. 

Unfortunately, with the absence of gravity, a negligible internal 
energy dissipation in the material, and no possible energy loss with the 
non atmospheric environment, there is no natural phenomenon that can 
provide damping for the vibrations and assure stability. 
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All these characteristics xmke active dynamic control of a flexible 
structure compulsory, so that vibrations and geometric distortions be 
kept within an acceptable level. The National Aeronautics and Space 
Administration (NASA) is scheduling a series of experiments, referred to 
as Control Of Flexible Structure, or OOFS program, in order to 
investigate, and to validate a technology data base for the suppression 
of inherent dynamic responses in large flexible spacecraft aind the 
avoidance of undesirable interaction between flexible structures and 
controls. The OOFS I experiment will be a 60 meter high truss mast 
bearing a 100 kg tip mass, and will be deployed from the cargo bay of 
the Space Shuttle. Figures 1-1 and 1-2 illustrate its characteristics in 
greater detail. The mast structure is made of more than 600 struts put 
together with about 1200 pinned-joints whose characteristics are 
described in figures 1-3 and 1-4. This study is a contribution to the 
OOFS project, and is aimed at obtaining a better understanding of the 
effects of the joints on the structure dynamics. 

1.2 Review of Previous Related Works. 


Studies have already been conducted concerning the characteristics 
of the joints and their effects on the structure dynamics from a linear 
point of view, assuming that each joint has a constant stiffness and a 
constant viscous daunping, or from a nonlinear point of view, assessing a 
more realistic behavior of the joints. 

Among recent works, R.Y.-K. Lee, [1], examined the linear effects 
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of the joints on a three piece boom. He also derived nonlinear models 
for the joints, based on a physical understanding of what causes the 
nonlinear behavior, such as Coulomb friction or surface mating. These 
models were then used to simulate dynamic responses. Problems in 
simulation raised by the nonlinear character of the systems were also 
investigated for various integration schemes. 

K.W. Belvin [2] also presented a thorough investigation on how to 
model the nonlinearities in the joints, how to obtain a simpler model 
more useful for larger problems, auid how to use these models, along with 
an elaborate finite element analysis of a truss structure allowing large 
deflections, to simulate its dynamics. Next, he analyzed some of the 
effects of the joints on the dynamics of a plane truss structure with 
the new refined simulation technique, and obtained some interesting 
qualitative indications of its behavior. He also addressed the problem 
of parameter identification, eind derived a method to fit optimally , sin a 
least-square error sense, his model to experimental test results. 

Parameter identification was the only goal of the thesis by 
K. J. O’Donnell , [3], where dynamic characteristics of the joints are 
found through force-state mapping techniques which yield very realistic 
and accurate models. 

1.3 Thesis Objectives and Organization. 

All of the aforementioned studies look only qualitatively by meeins 
of simulation at the effects of the joints on free structure dynamics. 
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They do not address the problem of the structure behavior under the 
effect of a control system, nor do they take into account the 
applicability of the simulation techniques to very large structures. 

This present study is therefore aimed at filling this gap and at 
investigating some possible approaches to analyzing the behavior of 
nonlinear truss structures under active control. 

Chapter 2 shows how to model the problem, and specifically how to 
treinsform the continuous mechanical problem into a finite state variable 
representation which is very well suited for the design and the study of 
control systems. 

It would have been possible to investigate numbers of th system 
properties using the newly developed finite state description, but among 
the different issues, the major problem of stability was not addressed 
and the study was focused instead on the possible occurrence of limit 
cycles in the structure under active control. In fact, stability appears 
to be easy to guarantee for the particular type of problem found with 
the control of the COFS I Mast, and it is felt that the design of a 
controller that stabilizes the structure should be a rather simple task 
in this case. On the other hand, significant nonlinearities can interact 
with the control system emd produce sustained oscillations, also called 
limit cycles, and because of the expected properties of the joints, 
those limit cycles seem likely to occur. The reasons for particularly 
studying limit cycles are developed more detail in Chapter 3, with the 
rest of the chapter being devoted to deriving analytical tools to 
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predict their occurrence. 

Limit cycles are studied under a restrictive single harmonic 
hypothesis. The h3^othesis allows the nonlinearities in the joints to be 
talcen into account through the use of describing functions, giving means 
to derive methods applicable to large systems for the search of limit 
cycles and for the eventual determination of their stability. The 
special form of the (DOFS I Mast problem seems to Indicate that the 
accuracy of the solutions found through these methods may not suffer too 
much from this simplifying hypothesis, and their use seems, therefore, 
to be rationalized. 

In (Dhapter 4, the practical aspects are taken into account, and the 
section shows how the theoretical methodology of the previous chapter 
can be computationally implemented to determine limit cycle existence in 
real problems. Minimization methods are used in the limit cycle search 
algorithms derived. Some of the presumedly most effective minimization 
methods are reviewed in the section, and their pros and cons are 
discussed. The computational task associated with the different 
numerical implementations is evaluated in order to facilitate the 
selection of a most efficient algorithm, and also to permit the overall 
evaluation of the effectiveness of analytical methods, and ultimately to 
allow to compare it to the effectiveness of alternate methods. 

Finally, some simple examples are derived as an illustration in 
(Dhapter 5. They allow comparisons between the different numerical 
techniques for very small order problems, and they confirm results 
previously obtained on their efficiency. 
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The conclusions that can be drawn from this first very general 
study about the applicability of analytical methods to predict limit 
cycles in very large structures under active control are presented in 
Chapter 6. 

Possible improvements of the analytical methods, as well as 
alternate methods which should constitute the object of further studies, 
are indicated in this last section. Among the alternate methods, the 
determination of limit cycles through simulation is more specifically 
detailed, and a rough estimate of the computational task associated with 
it is given. The comparison of the estimates of the calculation 
requirements obtained for the two different approaches, the analytical 
approach and simulation, gives a strong indication that simulation 
should be in fact one of the most effective approaches to the 
determination of limit cycles in large structures comparable to the CX)FS 
I Mast. It is therefore believed that simulation techniques for large 
structural systems should be further studied in order to improve them 
and to reduce the computation they require, which still appears to be 
considerable when dealing with large dynamical systems like space 


structures . 


r 

r 

r 


17 


origikae psge is 

OE POOR quality, 


m (QUAU^-Y 





18 


y (tip 

neter 
nsors , 







20 




21 


CHAPTER 2 


PR(»LJEX NCXKLING 


2.1 Basic Hypothesis on the Nature of Problans in Structural Mec^aTiir.w . 

In most of the developments of Structural Mechanics, the hypothesis 
is made of the linear behavior of continuous materials. The strain, or 
local distortion of a body, is linked to the stress, or effort, through 
linear relations which use characteristic parameters of the material 
like the Poisson coefficient and Youth’s modulus for Ein isotropic 
material. This h3rpothesis is verified as long as the body is not subject 
to large deformations, and allows a great simplification of the study of 
the djmamics of the non-rigid body. 

2.2 Review of a Conmon Approach for the Analysis of the Dvnamirs 
of Non - Rigid Bodies. 

2.2.1 Problem Statement: the Modal Analysis. 


The dynamics of a non-rigid body are governed by a set of partial 
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differential equations, along with a set of boundary conditions. Under 
the linear hypothesis, the general form of the problem can be written as 
follows ( Meirovitch [4]): 

Lu+Mu=f ,P€D (2-1) 

B. u = 0 . P € S: i=1.2....p (2-2) 

1 ~ ~ 

where L is a differential operator matrix of order 2p, M is a mass 
matrix, u(P,t) is the displacement vector at time t of the point P of 

i» 

the domain D occupied by the structure, u being its second partial 
derivative relative to time, B. *s differential operator matrices of 
strictly smaller order than L, and S the boundary of D. 

£(P,t) represents the external loads which can be either distributed or 
discrete, where in the latter case it can be expressed as: 

f(P,t) = F(t)6(P - P^). (2-3) 

where 6(P - P^) is a spatial Dirac delta function. 

Because of its linear character, the problem can be solved in a 
particular way, through modal analysis: periodic solutions of (2-1) are 
sought when no external load is applied. Exponentially decaying 
solutions will be sought if daimping is introduced, but the structure’s 
internal damping is usually assumed to be zero. A mode r is then 
represented by a complex pareimeter that describes the frequency of 
the solution and a mode shape (j) which is a function vector that 

/v#r 

associates the maximum amplitude of the displacement to each point of 
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the body. Looking for solutions of the form 

exp(X^.t) 


(2-4) 


in equations (2-1) and (2-2) yields the following equations for (<b ,X ): 

..Lr r-' 


(L + M Xj ) ^^ = 0 . P € D (2-5) 

®i $r = 2 *P ^ S .U1....P (2-6) 


which is an eigenproblem ( [4]). 

The eigenfunctions (|) ’s constitute an orthonormal basis of the 

space of the solutions, provided they have been normalized, and for a 
chosen dot product, they satisfy: 



for any two modes r tind s, where 6 is the Kroenecker symbol. 

rs 

Solutions for any kind of loading can therefore be projected on the 
space spanned by the modal solutions. Writing the solution u as 

CO 

u(P,t) = 2 u (t) 6 (P) , u (t) C R , (2-8) 

r=l ^ r 

transforms the partial differential problem into an infinite dimensional 
system of ordinary time differential equations satisfied by the modal 
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coordinates u (t): 

u - u = f (t) . r=l...«>. u € K (2-9) 

r r r r 

where : 

f (t) = \ dD (2-10) 

is the forcing term for mode r. 

2.2.2 Practical Solution for a Real Problem. 

To find the eigenvalues and eigenfunctions of the modes is usually 
very difficult and closed-form solutions only exist for very particular 
problems. Approximate solutions for a limited number of modes can 
however be found through finite element analysis for example. 
Approximate eigenfunctions (fP associated with approximate eigenvalues 

are defined by a finite number of parameters or degrees of freedom, and 
thus only s|>ELn a finite dimensional subspace of the solution space. The 
problem is discretized and can therefore be treated computationally, 
whereas the continuous problem cannot be treated at all. 

Discretization preserves the orthonormality property and the 4P’s 

still constitute an orthonormal family on which the general solution can 
be projected. 

The error of the approximate solution can be timed by choosing an 
adequate number of degrees of freedom to describe each mode shape. 
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However, the number of degrees of freedom has to be increased if a 
better accuracy, and a larger number of modes, are wanted. 

A practical solution will therefore be sought in the form ([4]): 


r=n 


u(P.t) = 2 u (t) (}f(P) 
r=l ^ 

where n is the number of selected modes. 


( 2 - 11 ) 


The modal coordinates u^’s satisfy the following finite dimensional 
differential system: 


with: 

2 ind: 


I . r=l...n, 

f 

= € (P)-f(P.t) dD, 

Jd ~r 


( 2 - 12 ) 


. a 

A = J CO 

r "" r 


= -1 


A lumped linear system can therefore approximate the dynamics of 
non-rigid bodies, simplifying tremendously the analysis, even though the 
dimension of the linear system must be very important for a satisfactory 
accuracy to be reached. Any solution to the dynamical problem for the 
body will be sought only in the form: 


u(P,t)=2'' u(t)4>(P) , 

r=l ^ 


where (A ,<t> ) is a mode or an approximation we have of it. 

I 


( 2 - 13 ) 
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2-3 Modeling of a Truss Structure with Nonlinear Joints. 

2.3.1 General Approach. 

Nonlinear systems with two distinct parts, one linear and one 
nonlinear, constitute a special class of nonlinear systems which have 
received a lot of attention already. Their general architecture is shown 
in figure 2-1 . 

Extensive study of the behavior of this type of system with a SISO 
nonlinearity is made in Gelb and Vander Velde [5]. Extension to the MIMO 
case exists in [6,7]. 

A truss structure is essentially constituted of linear elements, 
the struts, connected by nonlinear elements, the joints, and it would be 
appealing to put the problem under the aforementioned form to make best 
use of its linear characteristics. 

The most staightfoward approach would be to describe each strut 
with a linear model and cascade and feed these linear subsystems back 
through the nonlinear joints, but this would result in a maze of 
inextricable connections. Thus our approach is to first replace the real 
joints by linear equivalents. The resulting altered structure becomes 
linear, and modal analysis of the whole linearized truss is applied as 
described before. The linear model of the structure includes state 
variables for the joints which will be used to feedback the remaining 
nonlinear p>art as forcing terms. 
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2.3.2 Mathematical Formulation. 

Let 1,2,...,N. be the indices of the joints in the structure. 

The joints are assumed to be massless and each joint is described by 
only two nodes at its ends located at P.j^ and P ^2 • 

The efforts applied by joint i on the rest of the structure are 
considered discrete, and can be written: 

F.(P,t) = F^(P.t) + F^(P.t) , (2-14) 

where F^ is linear in the displacement of the joint and its rate of 
displacement and the remaining nonlinear part. 

The way the F^’s are chosen is theoretically unimportemt, as long as it 
yields a linear model for the entire structure: the nonlinear feedback 
will anyway cancel out any unrealistic dynamics assumed for the joints. 

It would be much more satisfactory, however, if the linear 
equivalent retained as much as possible of the real behavior for 
principally two reasons: 

First, it might allow to design the controller using only the linear 
part, when the remaining nonlinearities are small enough, and to apply 
therefore only linear system control design tools. 

Secondly, and more importantly, the modal analysis fixes the subspace 
where solutions will be sought. Thus, a main concern is that the 
component of the complete solution obtained on this subspace may not be 
preponderant. Perpendicular components may not be negligible, as it was 
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implicitly assumed, if the linear equivalents differ too much from the 
real joints . 

Throughout the study, we will assume that such a close linear 
description of the joint exists as suggested by the test curve we have. 
( figure 1-4). 

Since the linear part has been included to perform the modal 
analysis of the overall structure, the problem can be written as 
follows: 


" 2 

U + <0 u 
r r r 


2 J 

i=l 


+ f®f (t) 


t) * 

r=l,2 . . , 


n 


(2-15) 


6X t 

where f^ is an external forcing term, basically control forces and 
disturbances, and can be the result of either distributed or discrete 
loads. 

The matrix formulation of the problem is the following: 


+ A + B y + L.d 

~ J~NL c ~c ~ 


where ^ is the vector of modal coordinates 

s = C “r “2 ■ 


(2-16) 
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and 

2 2 2 

A = diag( Wj, ^2 ) . 

~NL vector of nonlinear forcing terms. 


1 1 J J 


(2-17) 


the B . matrix is: 
J 


"j = 


4>{(Pi ) - <l>{(Pi ) 

/VpA. ± ^ /v»A ^2 




Jl J 

J «] 


^Cn.!) - ^C’n. 2 ) 

3 J 


(2-18) 


a vector of supposedly discrete control forces applied at 

Q 

different locations P.’s, 

J 


B = 
c 


F^CP^.t) 

El(Pr') ■ 



$If2) 1 ■ 

■ if[‘0 

1 

gcp^) 

1 

i 

1 


1 

^(■■ 2 ) I- 

1 

• • i 4>^(p‘") 

1 Zn^ 


(2-19) 


( 2 - 20 ) 
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d represents the disturbances and L the matrix through which they act on 
the system. 

We implicitly took massless joints by assuming that the sum of the 

efforts on any joint is zero. Hence, the term contains only the 

forces applied at one end of each joint. 

The nonlinear terms depend on the joint’s state variables, and we 

can write as : 

~NL 

^isn = ^MT <S) 

~NL ~NL'-~ ~ z 

where q contains joints’ history states which are necessary to 


completely describe the joints’ behavior, and whose dynamics are usually 
nonl inear . 

The block representation of figure 2-2 summarizes the way the open- 
loop system is described. 


2.3.3 Dimensionality. 


All vectors and matrices presented before have various dimensions, 
and it is important to know them to appreciate the size of the problem 
and later discuss the computational load and the size of memory required 
for its resolution. 

The modal vector ^ is an n x 1 vector and the modal matrix A is 
therefore n x n . 

Each eigenfunction (j) is a 6 x 1 vector function, since it is 
necessary to include rotational deformation along with displacement in 
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the finite element representation. 

The discrete forcing terms F are also 6x1 and include a three 

dimensional force as well as a three dimensional torque. 

Hence, is a {6 N.) x 1 vector, and B. isan nx(6N.) matrix. 

~NL ^ J J 

Similarly, ? is a (6 m) x 1 vector, and B is an n x (6 m) 
'^C c 

matrix, where m represents the number of actuators. 

Typically, the OOFS I experimental mast has about 1200 joints, aind 

we will take thereafter N. = 1000 as the order of magnitude to discuss 

J 

the size issue. The modal analysis should therefore include at least 
1000 nodes and about 10,000 degrees of freedom. The number of modes to 
be retained should be in the range of 100. Hence, the model yielded by 
the modal analysis is considerably reduced compared to the one used to 
perform the analysis. However, it should still be meaningful if the 
modes are carefully selected. 

2.3.4 Introducing the Control Laws. 

As mentioned in the introduction, active control is required on 
most of the large flexible spacecraft to transform their dynamics and to 
obtain an acceptable behavior. 

The type of control laws we will consider here are linear control 
feedback laws. It might be, for example, full state feedback control, or 
a LQG/LTR controller: the main purpose of this research is not, however, 
to study the design of controllers for linear descriptions of large 
space structures, and we do not include in this discussion the use of 
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reduced order dynamic controllers nor a wide variety of design methods 
which appears in the literature. 

Thus, the control vector has the general form: 

= ~ ^1 31 - S ^ • (2-21) 

where and G 2 are two (6 m) x n gain matrices. We CEin therefore 
rewrite the problem as follows: 


J -f B_, j -f {A + Gj) J = Bj + L d . (2-22) 

The state vector can also be augmented if integrators are to be 
included in the control system, before the actuators’ command inputs for 
example, in order to have a zero steady-state error in the presence of 
constant disturbances. This operation is classic in linear feedback 
control design ( [8]), eind the resulting form of the problem is similar 
but with an increased state vector. 

As we can see in figure 2-3, the system still keeps a desired form 
with separate linear eind nonlinear parts when the control loop is 
closed. Closing the loop only affects the linear part, and there is no 
need to repeat the modeling process for different control feedbacks. 

In summary, only one very expensive modal analysis of a linear 
equivalent of the entire structure with no active control is required, 
yielding a linear model whose dynamics is modified by the introduction 
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of various feedbacks, one nonlinear corresponding to the remaining 
nonlinear characteristics of the joints, eventually a second nonlinear 
one corresponding to a possible nonlinear control feedback law, and a 
linear one corresponding to the linear control feedback law. As the 
control feedback loop is closed on the linear dynamics of the open-loop 
model, a new closed-loop linear system is formed, and this system 
becomes the linear part of the total system, to which the nonlinear part 
is fed back. The model of the system stays, therefore, in the form of a 
nonlinear feedback system, which allows simpler and further studies of 


its behavior. 
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CHAPTER THREE 

ANALYTICAL PREDICTION OF LIMIT CYCLE 
IN LARGE TRUSS STRUCTURES 

3.1 Foreseen Importance of Limit Cycle in Behavior of Large 
Truss Structure under Active Control. 


Nonlinear systems can behave in a lot of different ways, and any 
designer should look for properties such as the presence of multiple 
equilibrium points, their stability, the possibility of jump phenomena 
and the occurrence of limit cycles and amplitude dependent behaviors in 
the closed-loop system. Everything is equally important for the 
performance of the design, but some assumptions made about the joints in 
the OOFS I mast make some of them less critical. 

3.1.1 Multiple Equilibrium Points. 

Multiple equilibrium points means the possibility of having steady- 
state geometrical deformations of the structure and thus a loss of 
conformity in its shape. It could be very dsimaging if the geometry is a 
preponderant factor in the completion of the mission, as it is for large 
flexible antennas. 

Even though the study does not address the issue directly, we will 
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show how to modify the limit cycle search technique to be able to 
determine the equilibrium points. 

3.1.2 Stability about an Equilibrium Point. 

Stability is always a major issue in any control problem, and if it 
is quite straight foward to evaluate it for linear systems, special care 
has to be brought when it comes to nonlinear cases. The center meuiifold 
theory permits one to go further and extend linear system stability 
criteria to nonlinear systems to assess local stability, and Lyapunov’s 
second method is a very powerful tool used to answer global stability 
questions. ( [9] and [10]). 

However, in the case of the joint dominated truss structure some 
heuristic arguments will easily convince one that stability can be 
guaranteed provided the joints possess certain properties. More 
precisely, we suppose the joints to be asymptotically linear, or if f 
represents the load at one end of a joint and x its displacement, we 
suppose that there exists a stiffness K such that: 

lim (f(x) - K X )/f(x) = 0 . (3-1) 

X -» « 

The range of x is naturally limited inside the region where the 
materials have a linear behavior, but it can be assumed that this region 
is wide enough, compared to the area where the nonlinear effects are 
predominant, so that (3-1) is satisfied. The property is illustrated in 
figure 3-1 . 
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Another underlying assumption made here is that the joints have no 
real d 3 mamical characteristics, or that the restoring force depends 
mostly on the displacement x euid not the rate of displacement x, and 
that the energy dissipated per cycle depends on the eunplitude of the 
excitation rather than its frequency. Otherwise, a viscous damping term 
must be added to the asymptotical linear model in (3-1). 

The pinned- joints used in the OOFS experiment seem to fall in this 
category of asymptotically linearly behaving elements because the main 
nonlinear phenomenon is thought to be the play at the pin connection suid 
the Coulomb friction which opposes the backlash. 

If the linear asymptotes are used to build the linear equivalent 
model of the structure, and a controller that stabilizes it is added, 
then the closed-loop nonlinear system will not blow up, since for large 
enough initial conditions the joints will behave like their asymptotical 
equivalents and the state variables will decay. Hence, bounded initial 
conditions result at least in bounded responses, which means that the 
system is stable. 

3.1.3 Jump Phenomena. 

Jump phenomena occur when the system switches from one equilibrium 
point to another for an infinitesimal change in one partuneter. They are 
also referred to as catastrophe and are mostly studied in the very 
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powerful theory known as the catastrophe, or bifurcation theory ([9])- 
Buckling is an example of it, and it is the usual jump phenomenon found 
in structural problems. 

Such things do not seem likely to happen in our system since the 
purpose of the structure is not to carry axial loads, and the presumed 
asymptotically linear behavior of the joints should guarantee continuity 
of the response of the structure. Hence, further study seems not 
important to this point. 

3.1.4 Limit Cycles and Amplitude Dependent Behavior. 

3. 1.4.1 Limit Cycle Definition. 

A limit cycle is a sustained, or self excited oscillation in a 
nonlinear system. Its representation in state space is an isolated 
closed path , and if the initial conditions are set on one point of the 
path, the trajectory of the system will remain on it. 

A limit cycle is said to be stable if the trajectories originating 
in its vicinity tend toward it as time goes on. On the other hand, it is 
said to be unstable if the trajectories go away from it. Figures 3-2a 
and 3-2b illustrate the concept of stable and unstable limit cycles. 

Unstable limit cycles are also called stability boundaries, limit 
cycle implying stability in that case. 

Self excited oscillations might occur in the structure because of 
the nonlinear character of the joints: the nonlinearities may introduce 
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some delay in the transmission of the actions of the controller and they 
may change the phasing to such a point that the controller would add 
energy instead of dissipating it. 

3. 1.4.2 Amplitude Dependent Behavior. 

Amplitude dependent behavior happens in the presence of two limit 
cycles and is best explained in figure 3-3: 

The inner limit cycle is unstable, or is a so-called stability boundary, 
and the outer limit cycle is stable. Therefore, for initial conditions 
within the stability boundary, the system will come back to zero, and 
will be stable, whereas it will be trapped and will display a limit 
cycle if the initial conditions are anywhere outside this boundary. 
Hence, the behavior of the system depends on the amplitude of the 
initial conditions. 

3. 1.4.3 Consequences of Limit Cycle Occurrence. 

Limit cycles in a structure can cause trouble of the same nature as 
free oscillations do when no active control is applied. 

Their worst effect would be to interact with dynamics of modules or 
other systems connected to the truss and excite those dynamical systems 
periodically near their own resonant frequencies, thus forcing an 
amplified response or even causing instability. 
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But even if unstability is not at staJce, vibrations can be 
detrimental for the completion of the mission: loss of accuracy in the 
geometrical shape can result in a less efficient reception of solar 
power for solar arrays, or it can degrade the signals received and 
transmitted by flexible antennas. The vibrations can induce an 
unacceptable level of acceleration for microgravity experiments taking 
place in a lab facility of the Sp>ace Station. 

Finally, it can be very damaging for the spacecraft itself, and it 
can reduce its reliability and its mission time by speeding the process 
of fatigue in the joints, in the electronics on-board and usually very 
sensitive to vibration, and also in the sensors and actuators of the 
control system which interact with the structure to produce the limit 
cycles and which are continuously sollicitated. 

Therefore, the control system should be designed in such a way that 
no possible oscillation can occur, and searching for limit cycles must, 
hence, be of prime concern. 

3,2 Analytical Determination of Limit Cycle: Problem Formulation. 

3.2.1 General Approach, 

A.A. Aderibigbe made a thorough study on how to predict limit cycles 
in multivariable nonlinear dynamic systems ( [11])* He reviewed some 
already used techniques and he presented an alternate one based on a 
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particular Ritz-Galerkin method known as the Harmonic Balance Method. 

From his work it appears that all applicable analytical techniques 
share some common f eatures : 

First they reduce the general multivariable nonlinear problem using 
quasi-linearization techniques in order to get a practical model to work 
with. As we will see later, quasi-linearization methods involve writing 
the dynamic equations as linear equations with coefficients that depend 
on certain properties of the state variables and which are chosen in 
order to reduce in some way the error made in the approximation. 

Then a search technique is employed. It is usually an algorithm 
which will iterate until limit cycle conditions are reached. It also 
usually tries to make best use of the quasi-linearized form of the 
problem and track eigenvalues or characteristic equation roots. 

The effectiveness of a technique can be appreciated by the degree 
of approximation reached in the process of getting a workable model and 
by the efficiency of the search algorithm in terms of memory and time. 

3.2.2 Quasi-Linearization Methods. 

Since a limit cycling behavior is essentially oscillatory, 
periodicity in the dynamical equations has to be considered; the quasi- 
linearization processes perform a sort of Fourier series expeuision. 
However, there are different approaches to do it, involving global or 
local approximations and retaining only one or more harmonics in the 
model. Three basic methods can be used. 
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The first method uses a global statistical linearization. Assuming 
the system is in the form: 

x=f(x)+Fw (3“2) 


where x is the state vector which contains a bias x, , a sinusoidal 

ru> ~'D 

component x^, and a raindom component x^: 


X = X, + X + X 

~ 


(3-3) 


£( x) represents the system’s nonlinear dynamics, and w is the input 
with: 


w = w, + w 

/v» 'V* 


(3-4) 


where w, is a bias input and w a random one. 

f^( x) can be approximated by the quasi-1 inear form 

f( x) = N.x, + N X + N X 

s^s r~r 


(3-5) 


,N^ and are obtained by minimizii^ the mean square error of the 

difference e 

e = f( x) - - N X - N X 

S’^s r^'^r 


(3-6) 
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The general form of the result is shown in Spanos eind Iwan 
( [12]). and is the following: 

(3-7) 
(3-8) 
(3-9) 


Nb - (E[f sj]) 

= (E[f sl]) 

N, = (E[f /]) {E[x^Js^])-‘ 


An alternate method consists in replacing individual constituting 

elements by their own quasi-linear representations. If x(t) is the 

scalar input to the nonlinearity, it is assumed that it only contains a 

bias X, , a sinusoidal term x , and a random one x : 

D s r 

^ ^ ""r 

The output is also assumed to be the sum of a bias, a sinusoid and a 
reindom signal and is written as: 

y(t) = + n^x^ + n^x^ (3-11) 

The gains and n^ also minimize the mean square error of the 

difference between the real and the approximate signal and are called 
describing functions ( [5]). It is shown in [12] that for systems 
consisting of isolated nonlinear elements connected between nodal 
points, the global and the element-by-element quasi- linearization yield 
the same model. 

A third method uses the harmonic balance method and describes the 
nonlinear system as a proportional plus derivative system as a whole. 
Its main feature is that it involves more than one harmonic in the 
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assumed form of the state variables, as detailed in [11]. 
The assumed form of the state vector is: 



where the subscript b stands for bias, and the subscript s^ for 
harmonic of the sinusoidal term. 

The system : 

tt 

X = f(x , X ) 


(3-12) 

the 


is rewritten using: 


ffx .X ) = N,x, + N -X + - Nj.x + 
^ / b^ pl^s^ (i) dl~s^ 


1 


'1 


+ N V X +, - Nj, x 

ph'^Sh h(j dh~Sj^ 


(3-13) 


where each harmonic term and are again found to minimize a mean 
square error. For only one single harmonic, the solution yielded is the 
same for all three approaches. 


3.2.3 The Single Harmonic Hypothesis. 

In simulations performed in [2], a four bay plane truss structure 
excited by a sinusoidal input responded with a marked sinusoidal 
character and without any other apparent harmonic terms. In the presence 
of active control, it is also presumed that higher harmonic resonsinces 
will be daimped enough so as not to appear in the response of the 
structure. We must emphasize here that the linear system to which the 
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nonlinear joint effects are applied is the closed-loop controlled 
structure. Therefore, the structural resonances within the bandwidth of 
the control system will be damped, and the bandwidth of the overall 
system will be limited. Hence, it is believed that a study of possible 
sustained oscillations involving only the first harmonic will be enough 
to predict whether or not limit cycles can occur. Anyhow, an analytical 
study involving higher harmonics would result in the multiplication of 
state variables and would rapidly make the problem untractable. 
Furthermore, no stability results can be derived if more than one 
harmonic is retained. 

We shall mention however that an attempt to quantitatively evaluate 
the quality of the first harmonic approximation is presented by 
Qiouldury and Atherton ( [13]). It involves the computation of the 

distortion of the signal fed back after a sinusoidal input is sent in 
the open-loop. If the signal x(t) has a Fourier series expansion of the 
form 

00 

x(t) =2 a,cos(kut) (3-14) 

k = 1 

the distortion is 

« 2 1/2 

A = '■k=2 ^ . (3-15) 

and it is claimed that if the distortion is less than 6%, one can expect 
an error no more than that amount. Hence there exists at least one test 
that can be performed to validate the approximation, but it requires 
mostly simulation techniques to be performed. 
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3.2.4 Quasi-Linearization of the Problem using Dual Input Describing 
Fimctions. 

Because of the form of the problem, quasi- linearizing the nonlinear 
part of each joint is obviously the best way to operate. Dual Input 
Describing Functions (DIDF) are used under the single harmonic 
hypothesis. 

3.2.4. 1 Dual Input Describing Functions. 

DIDF involves SISO nonlinearities. If x is the input and y the 
output of such an element, x is first assumed to be: 

X = B + A sincjt (3-16) 

where A an B are real, and y is related to x in the most general case as 
follows: 


y = ng(A,B,co) B + n^(A,B,(j) A sin<ot , (3“17) 

where n^ is a real, and n^ is a complex, gain chosen so that the mean- 
square error of the difference between the true and the assumed signal 


is minimized. 
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As shown in [5], the general forms for the gains are: 


1 


= I y{ B + AsinG, AgjcosG )d0 

^ 2ir B *^0 


(3-18) 


"a = ^ J ■ 


n 

P TT A ^0 


y( B + AsinG, AucosG ) sinG dG 


(3-19) 


1 

n = y( B + AsinG, AucosG ) cosG dG 

^ TT A 'Jq 


(3-20) 


A simplification occurs when it is supposed that the output y only 
depends on the input and not its first time derivative. If 

y = y( X, x) = y( x) (3-21) 

then there is no term in AucosG in the expressions 3-18 to 3-20. The 
same simplification occurs when, as it is assvuned for the joints of the 
CX)FS I Mast problem, the nonlinear function is symmetrical, even though 
not necessarily single-valued. 


3. 2. 4. 2 General Quasi-Linearized Form. 


The inputs to the global nonlinear part of the system are the modal 
coordinates u^’s contained in the modal vector V. The outputs are the 
forcing terms applied at one end of each joint, i.e. the 's. The 
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global representation of the nonlinear part has the following form: 


^NL = ^ h' 2 * 1 ) exp{j(ot) (3-22) 


where the modal vector is supposed to have the form 


where : 


S = ^ + ^1 exp(jo)t) 

^ is a n X 1 real vector representing bias 


(3-23) 


is a n X 1 complex vector representing the amplitudes as 
well as the phases of the oscillations of the state variables. 




^ and are both (®^j) ^ ^ matrix functions. 


3-2. 4. 3 Quasi-Linearization Process for One Single Joint. 


The first step in getting the global model is to find the quasi- 
linearized model that gives for one joint i the forces in terms of 

and . 

Mj ~1 

The use of DIDF implies a SISO relation between an excitation term 
and an output from the nonlinear device. However a joint is described by 
6 parameters: 3 displacement parameters and 3 rotational parameters, aind 
the efforts it transmits are also 6 dimensional, consisting of a 3 
dimensional force eind a 3 dimensional torque. This indicates that 36 
coefficients should be found to fully describe the behavior of the joint 
in the most general case. 
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Hopefully however, there should not be too many nonlinear effects 
and they should be decomposed into simple SISO nonlinear relations such 
as the one relating the axial load to the axial displacement in figure 
1-4. Hence, the decomposition of the problem will have the following 
steps. 

First, consider the vector of parameters called that defines the 

deformation of the joint: 

i T 

q=[vvv000] f3-24l 

^'■xyzxyz-' 

^x’ ^y ^2 components of the deformation of the joint 

measured at its end 3.1ong the reference coordinate frame axes Oxyz, 

®x* ®y ®z rotational distortion of the joint’s 

body. The deformation components in a two dimensional case are shown in 
figure 3-4. 



Figure 3-4: .Joint’s Deformation Parameters 
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is obtained by a linear transformation of the state vector 

/V/ 

and it can be written as* 


= C. 

where is the following matrix: 


( 3 - 25 ) 






( 3 - 26 ) 


Now suppose that one particular nonlinear effect has been retained, 
which depends on a one dimensional quaintity which is linearly related to 
q^ , and whose output is also one dimensional. The direction along which 

the input and the output respectively occur are P? and P^. The input 
to the nonlinear element is therefore: 



( 3 - 27 ) 


X can be separated as the sum of a bias B and a sinusoidal term 
of amplitude A and phase which are directly related to the bias 
and sinusoidal components assumed for the modal variables as shown: 


= £i ( 3io ) 

T T 

^ ^1 (3-28) 


= B + A exp( j(ut + ♦)) 


The one dimensional output y is approximated using DIDF by: 


y = ng(A.B) B + n^(A.B) A exp( j(ut + (())) (3-29) 


= VI Cl Si 

+ nJ I P. C. % 
A' ' 1 ~1 


T T 

C. ^ c. <w_ + 


T T 

Ei S ^ ^ Ei S ~i exp(jwt) 


(3-30) 


where | z | represents the magnitude of the complex number z. The phase 
term (|) does not appear explicitly in (3-30) since is already 
expressed as a complex vector, and the sinusoidal term in x is therefore 


described with a complex sunplitude. 
The resulting forcing term is: 
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Hence, for one given SISO effect recognized in the joint, the load- 
displacement law can be written as follows: 


2 !i) So * "a< Si exp(j..t) (3-32) 


where ^ 


and 


T T T 

Ni = P® n„(l C. 'W. 1 , P^ C. ^ ) P^ C, 


n t'T 

N* = PT n.{ \ P: C. <U, \ . p! C. JIL ) P. C. 
A ~i A'^ ‘ ~i 1 ~1 ' ~i 1 ~i 1 


( 3 - 33 ) 


( 3 - 34 ) 


Nj and are two 6 x n matrices that depend on the characteristics of 
A B 

the joint through the direction of the main exciting variable, the 
direction of its action, the nature of the SISO nonlinearity summarized 
by a describing function and the way the exciting deformation of the 
joint is related to the state variables. 

It is possible to consider that more than one nonlinear effect 
exists, but in that case the process shall be repeated and the matrices 
as well as Ng found for each process summed up, provided all the 
nonlinear effects are SISO with a good approximation. 


3. 2. 4. 4 Quasi -Linearization Process for the Whole Structure. 

The general model is obtained by assembling the different models 
found for each joint. Since the general forcing term has the form 
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shown in (2-17), it can be rewritten in the form (3-22) with: 


VSo-- Si) = 




"b(3!o- ?i) 

"H- Si) 

"2'(So’Si) 

"I'So- Si> 
4%. Si> 


3.3 Limit Cycle Conditions. 


3.3.1 Closed-Form Conditions. 


(3-35) 


(3-36) 


A limit cycle occurs when an oscillatory regime is established in 
the system eind when both the nonlinear forcing terms tind the modal 
variables are sinusoidal functions of time. 

In steady-state, the state vector ^ is related to the forcing term 
~NL a linear matrix transfer function we call G, which can be 
derived from equation (2-22). Equation (3-22) on the other hand shows 
the relation between the forcing term and the modal vector <?/. 
Combining the two equations, and balancing the harmonic terms, leads to 
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a system of 2 n equations, with n of them on K, when the bias terms are 
considered, and n of them on (D when the sinusoidal terms are treated: 

[ I - G( 0) = 9 

- a ^ % e (3-37) 

~0 '^1 

[I - G(j(o) %Q. %^) ] = 0 

where I is the n dimensional identity matrix and where 0 is the n 

/v» 

dimensional null vector. 

The unknowns in (3-37) are ^ and the frequency (j. If a limit 

cycle exists, the system must admit a solution with different from 0, 
which meeins that there must be at least one modal coordinate that 
displays a sinusoidal oscillation of non zero amplitude. 


It must be mentioned here that an alternate system can result from 
the combination of equation (2-22) and equation (3-22) which results 
from the elimination of the state variables between the two expressions: 

[ I - So’ <=( 0) ] ?o = 9 

(3-38) 

. c I - V So- Si> G(J«) ] = 0 

where I is the N. dimensional identity matrix and 0 is the N. 

J ~ J 

dimensional null vector. 

Since it is supposed that N. is much greater than n, it is obvious 
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that the system (3-38) must contain redundant equations, and in einy 
case, the system of the smallest dimension must be retained. 

3.3.2 Expansion of the Equations on R. 


HEuidling complex variables is not a natural operation for 
computers. Therefore, in view of a numerical resolution of the system 
(3—37), and in order to be able to evaluate the computing task, it is 
relevant at this point to develop the equations on R. This requires 
first to express the variables on R as 


2o= C 


r ^2- 


0,T 

^ ] 
n-* 


^ 


1’ 



1 

u 

n 



(3-39) 


where the u^’s are all real variables and where the (|)j’s represent the 
phases between the sinusoids coming into each mode. The phases of the 
different signals are defined to an arbitrary constant, and the phase of 
the first sinusoidal component has been tedcen as zero without loss of 
generality. In a vector form, we can regroup the u^'s in 


Ho 

/ViX 


r 0 0 

[ u. , u„. 


r Ir Ir 
[ u. , u, 


2 ’ 


0,T 

• • "n^ 

lr,T 

li^T 

• • . Un ] 


Uj = [ 0 , u“ 


(3-40) 
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where f or any i : 


1 


u . e 
1 


Ir . . li 

= Ui + J u. 


(3-41) 


The transfer matrix G as well as the DIDF matrix shall also be 
sejjarated into their real Eind imaginary p>arts. Jil^ can be written as 


(3-42) 

A p q 

where jV and ^ are functions of and the frequency w. A 

p q MJ 

simplification appears when partitioning G that can also save an 
expensive matrix inversion if we recall (2-22). Expecting a sinusoidal 
solution in (2-22), the equation can be rewritten as: 

(-(j2 I + A + G^ + jo) G 2 ) (3-43) 

Therefore, (3-37) can be expressed as the following system of three 
subsystems of n simultaneous nonlinear equations 


[A t - BjAb(Uo,uJ.u})]Uo = 0 

[A-c^ub^^Gj - BjAfp(Uo.Uj.uj)]uJ- - B.J(^(Uo,uJ,uj)]u; =0 

[A-„2i^B^Gj - B^Afp(U^,uJ,uJ)]U;* [4,B^G2 - ,u5)]u;= 0 


(3-44a) 

(3-44b) 

(3-44c) 


The system can be solved in Uq, U^, and w, knowing that uj^ is 

null, or it can be solved in IL, the u^’s and the phase variables, but 

1 

the latter case requires the use of sines and cosines functions to 
further express the real and imaginary parts of the complex amplitude 
vector 
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3.4 Limit Cycle Stability. 

3.4.1 Introduction. 

Stability of limit cycles is an important issue since it can 
determine important stability properties for the system as seen in 
section 3.1.4. Simulation techniques are mostly used to assess stability 
by perturbing the system around a limit cycle. However an analytical 
determination can be made under the single harmonic hypothesis( [llj). 

The trajectory of a limit cycle is approximately an ellipse in the 
n dimensional state space, and may not be centered at 0 (figure 3-5 ): 

S - So - K'( Si eV"') 

% ~ %Q = Re( 2l2)cosut - Im( ^j^)siruiJt (3-45) 

In the subspace orthogonal to the plane IT^ of the ellipse, the 
system is believed to be asymptotically stable . The assumption is 
similar to the one used in the Hopf bifurcation theory ( [9]). and it 
states that the tangent linear model found about should have all but 
two eigenvalues strictly stable and only two on the ime^inary axis for 
formerly reviewed nonlinear phenomena such as a limit cycle to exist. 
The study of the nonlinear part is then carried in the 2 dimensional 
manifold e tangent to U at the limit cycle conditions. Therefore, the 
stability of the system in the subspace orthogonal to this plane does 
not affect the question of the limit cycle stability, and perturbation 
of the limit cycle for the purpose of testing stability should occur in 
to be consistent. 
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3.4.2 Consistent Perturbation Conditions. 


The parameters of system (3-37) are perturbed around the limit 
cycle conditions in terms of amplitude, phase and frequency. More 
precisely, if and define the limit cycle conditions, the 

perturbed terms will be 

J!o = So° 


% = + 6 % 

Ic 

X = jo) +6(7+ j6(o 


(3-46) 


If we rewrite the system (3-37) under the generic form 


^0^ 2!o’ ~i ’ ~ 9 


^x( %Q’ %i‘ - 9 


(3-47) 


the consistency conditions become 



(3-48) 


Furthermore, if the vectors of pareuneters are expressed in terms of 
their amplitudes and phases as in (3-39), then the condition (3-48) can 
be developed by linearizing (3-47) about the limit cycle condition. 
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n 6f 


6fo= 2 


i=l 6u 


+ 2 


^ 0 
6u. 
1 


Ic 


n 6f 


0 


i=2 64). 


64), 


Ic 


n 6f, 


+ 2 


i=l 6ut 
1 


6f 


0 


6(0 


6u!^ + 

1 

Ic 

(6u - j6ct) 


Ic 


(3-49a) 


n 6f, 


6fj = 2 


+ 2 


0 


i=l 6u. 

1 

n 6f, 


6u? + 


Ic 


i=2 64), 


Ic 


n 6f - 


i=l 6u. 

1 

6f. 


64). + — 


6(j 


6u. + 

1 

Ic 

(6(i) - J6 ct) 


Ic 


(3-49b) 


Regrouping the equations (3-49a) and (3-49b) in a matrix form, the 
condition (3-48) cein be written as: 


M 6V = -P 6c7 

/%/ 


where : 


6V = [6u^,6u2 6u^.6uJ.6u2. . 

. . . ,64)^.6<.f 


. ..60^.64)2,64)3. 


(3-50) 


(3-51) 


and where the matrix M and the vector P depend on the parameters at the 
limit cycle conditions. 

Therefore ([11]). a consistent perturbation can be put in the form: 


5V = -M ^P 6a 


(3-52) 



65 


3.4.3 Stability Conditions. 


The method to determine the stability of a limit cycle involves 
first some geometric trains format ions: basically one which trainsforms the 
coordinates so that the equation of the plane IT^ of the ellipse be in a 
canonical form, and that in the new coordinate fraune n— 2 states remain 
zero for any point of A scaling of the coordinates then brings the 
ellipse back to a circle. Details can be found in [11]. 

Stability cam then be assessed by applying a consistent 
perturbation which is also put in the new coordinate frame. Figure 3-6 
gives am illustration of the transformation and shows how stability cam 
be checked after that, as it is explained in the following: 

The distance to the origin of the point of initial conditions defined by 
the perturbation cam be measured to see if it is inside or outside the 
circle representing the limit cycle. The position is next compared to 
the sign of the parameter 6o which governs the size of the perturbation 
through the consistency condition (3-52): a stable limit cycle will 
require a positive 5a to place the perturbation inside the circle, and a 
negative one outside, whereas the correspondency is inverted for am 
unstable limit cycle. In other words, the sign of 


6a _ 5a 
5r ~ r - r 


(3-53) 


o 

is determined, where r is the distance of the perturbed point to the 
origin and r^ the radius of the tramsformed limit cycle. The sign must 
be negative for a limit cycle to be stable. 
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CHAPTER POOR 

NUMERICAL EETERMINATIOf OF LIMIT CYCLES 

4.1 Choice of a Resolution Method. 

4.1.1 Possible Appro£u:hes . 

The system (3-37) is a set of 3 x n simultaneous nonlinear 
equations. It is generally advised to approach the resolution of such 
systems by keeping their multidimensional vectorial form. using 
techniques based on the Jacobian of the sytem of equations. Numerous 
algorithms have been developed, the best of which being reported to be 
one based on a hybrid Newton-Raphson method [14]. A ready-to-use code, 
MINPACK-1 [15], is available and applies this method. The code is 
limited to problems with less than 20 variables. 

Dealing with the Jacobian becomes very expensive beyond that range 
of variables, and since the present problem is very similar to a 
singular value decomposition problem, it can be soved at least as 
efficiently with a function minimization approach. The methodology used 
here is therefore to try to minimize the square of the norm of the 
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residual vector e, where 


fn( Y ) 


fi( Y ) 


(4-1) 


and where f^, and V are the notations used in 3-46 through 3-51. 


4.1.2 Residual Functions. 


A limit cycle condition is reached if eind only if the set of 
variables satisfying (3-37), or its equivalent form (3-44), is such that 
is not the zero vector. Hence, using the variables as defined in 
chapter 3, and using the system of equations (3-44) rewritten in closed- 
form, the following residual function can be derived as 


R{ Y ) = 


II mo) - B .^g( V )]3 <q ir 

II % 11^ 

~1 

II cnco) - B V ir 

II % 11^ 

~1 


(4-2) 


where the matrix function is defined as: 


= -co^I + j(o B G„ + A + B G. 

C ^ C X 


(4-3) 


and where the rest of the parauneters are defined in chapter 2. Dividing 
the square of the norm of the residual vector by the square of the norm 


69 


of the complex amplitude vector ensures that any solution will have at 
least one non zero sinusoidal signal. 

Similarly, a residual function can be built to search for multiple 
equilibrium points. Since the origin is an equilibrium point, the 
residual function is straightfoward and is 

Ro( 3!o ) = " t So- 2 )]So (4-4) 

A set of limit cycle conditions equivalently satisfies (3-37) or 
sets the residual function R to zero. A limit cycle will therefore occur 
whenever the function reaches zero, and the parameters of the limit 
cycle are the ones for which this value is taken. Since R is always 
positive or null, it is also theoretically possible to infer the 
nonexistence of limit cycles by checking the absolute minimum of the 
function over the whole range of values the variables can take. If the 
absolute minimum is strictly greater than zero, no limit cycle can 
occur. On the other hand, if it is zero, limit cycles will occur 
whenever the residual function reaches its absolute minimum. 

It is clear, however, that restrictions apply, and a practical 
solution can only rely effectively on findlr^ a limited number of limit 
cycles, if some exist, since the computational task is the major 
limiting factor in the application of the theory. 
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4.2 Adapted Minimization Algorithms. 

4.2.1 Presentation of Ihicons trained Minimization Techniques. 

Unconstrained minimization can be achieved with various degrees of 
complexity which goes along with various effectiveness. It is obvious 
that the more information is known on the function behavior -i.e. the 
more derivatives we consider- the faster the convergence can be. 

However, the derivatives may not exist, but if they do, their 
computation can be very time consuming, and can require a lot of extra 
memory. Not only shall the methods with the highest rate of convergence 
be therefore considered to solve problems with a great number of 
variables. 

Almost all the minimization algorithms share the following 
features: from a current point, the algorithm finds a direction in the 
state space where the function is most likely to decrease. An univariate 
minimization, or line search, is then performed along this direction to 
locate the local minimum which becomes the new current point for the 
next iteration. 

The way those operations are performed separates the techniques, 
and it is possible to classify them into three main types: 

The first type is constituted by nongradient techniques which do not 
require any derivatives, but only grope toward the solution by 
evaluating the function at many locations, and make use of these 
computed values to determine the search directions. Their only advantage 
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is that they save the computation of the derivatives, but they suffer a 
very low rate of convergence and require many function evaluations. 

The second type is constituted by gradient methods which make use of the 
first derivatives to create a search pattern which yields search 
directions better distributed to find the solution more rapidly. Their 
rate of convergence is higher, but the line search is still blind. 
Finally, Newton and Quasi— Newton methods go one step further, and by 
using the Hessian or by approximating it, they allow the forecast of the 
location of the local minimum along a line search direction. Their rate 
of convergence is quadratic, but their main drawback is the necessity to 
update the Hessian matrix which can be a very large matrix. 

4.2.2 Possible Algorithms for the Determination of Limit Cycles. 

According to Scales [16] who reviews a great number of existing 
unconstrained minimization techniques along with their advanteiges and 
their disadvtintages, the Conjugate Gradient method is the most 

appropriate when the number of variables becomes very large, above 250 
variables. The main reason is that the method requires much less 
function evaluations than nongreulient methods, even if the gradient is 
derived numerically, but it does not require the computation of the 
Hessiain matrix which is square with dimension equal to the number of 
variables. 

However it is also reported in Scales that Quasi - Newton methods are 
very superior in terms of execution time and convergence rate. Among 
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those methods, the one suggested independently by Broyden, Fletcher, 
Golfarb and Shanno ([17,18,19,20]), and also referred to as the BFGS 
method is more specifically reported to be the most effective. The claim 
is also confirmed by Himmelblau [21], who conducted a uniform evaluation 
of the performance of different types of unconstrained optimization 
techniques. Even though the size needed to store the Hessian is a major 
problem when the function depends on a great number of variables, an 
other positive point is that the BFGS method suffers rather rough line 
searches, thus saving many function evaluations which appear to be one 
of the computationally most important burdens, as shown later on. 
Therefore, the BFGS method is another candidate to consider for the 
determination of limit cycles. 

4.2.3 Polak Ribiere Conjugate Gradient Method. 

The main feature of a conjugate gradient method is that the search 
lines are conjugate . For a non-degenerate quadratic function, this 
property means that the search directions form an orthogonal base, the 
dot product taken being the one induced by the Hessian of the function 
which is a constaint positive definite symmetric matrix in that case. If 

f( X ) = G X .X € R ™ (4-5) 

is the generic function to minimize, where G is a positive definite 
symmetric matrix eind m an integer representing the number of variables, 
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then the search directions p.’s satisfy: 


T T 

p G p = p G p 


5. . 
IJ 


(4-6) 


where 5. . is the Kroenecker symbol. Furthermore in the quadratic case 
the minimum of the function is reached in exactly m iterations. Things 
change in the nonlinear case since the Hessian is not constant. If m 
still represents the number of variables, the algorithm will perform 
major iterations, which are groups of m minor iterations,^ or 
computations of a search direction and an univariate minimization 


along the line. During a major iteration the directions are generated so 
that they satisfy the conjugacy property. It is proven, [16], that if g 

represents the gradient of the function f at a current point the 

direction Pj^ along which to find the next 


p>oint should be of the 


form 


Ek =-Sk^PkEk-i (4-T) 

to ensure conjugacy between the search directions in a major iteration. 

is a scalar which c£in be chosen according to the particular form of 
conjugate gradient method wished to be used ( [16]). The choice of 
due to Polak and Ribiere [22] reportedly succeeds more often aind should 
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be therefore preferred to other choices. Its form is 






T 



Sk-1 Sk-1 


( 4 - 8 ) 


The update of the search direction ^ requires the knowledge of the 
previous search direction the gradient at the previous iteration 

as well as the gradient at the current point 

When m minor iterations have been performed, the search direction 
is reset equal to the gradient at the current point: this reset has been 
reportedly very successful to prevent the search directions from 
drifting and lie in a subspace of k"* because of rounding errors in the 
update formula. 

Hence, only four m x 1 vectors have to be stored to update the 
search direction and this constitutes the main advsintage of the method 
for large systems. Also, when close to the minimum, the Hessian matrix 
is approximately constant tind a quadratic termination can be expected. 

The gradient can be found by finite difference as well as 
analytically without influencing the convergence properties of the 
method. However the method requires a very precise univariate search, 
tind each updated point coming from each minor iteration has to be very 
close to the true local minimum in order to get an acceptable 
convergence speed. This can be one of the major problems in applying the 
conjugate gradient method since the evaluation of the function can be 
very time consuming, and that is why other methods might be considered. 
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4.2.4 The BPGS Quasi-Newton Method. 

The Newton method uses Taylor series expeinsion of the function up 
to the quadratic term to find the next step, or equivalently a Taylor 
series expansion of the gradient to the first order. If the next point 
is related to the current point ^ as follows 

Xk+i = Xk * Ek (4-9) 

then the Taylor series expansion of the gradient g( X ) is: 

S' Xk * Ek> = «( Xk) * <=k Ek (4-10) 


where is the Hessiein matrix of f at A successful step is ope 

which sets the next gradient g( Xj^^j) to zero. Using the first order 

approximation for the gradient, this step is found to be 

Ek = - <=k‘*( Xk> • (4-11) 


which gives not only the direction, but also the step to take in that 
direction, at each iteration. 

The pure Newton method however suffers from multiple problems among 
which is the necessity to compute the Hessian G at each iteration and to 
invert it. which is expensive and not always possible. Quasi-Newton 
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methods have therefore been developed: they avoid the computation of the 
Hessian, and they carry and update instead an approximation of its 
inverse. 

As seen in [16], there exists families of Quasi-Newton methods, but 

the most successful seems to be the one referred to as the BFGS method. 

If represents the approximate of the inverse of the Hessian at the 
th 

iteration, it is updated as follows: 


«k*i = [I 


^ u [I - ^ ^ ^ 


(4-12) 


is the step taken in the iteration, the vector difference 

between the gradient at point ^ and ^ dimensional identity 

matrix and Hj^ the current approximation of the inverted Hessian. 

The step p, is given by 


where 


AX, 


Sk ‘ ^ 


(4-13) 


represents the step in a real Newton method. However the Quasi-Newton 

method starts only with an approximate Hessian, and therefore the method 
only considers Pj^ as a search direction and still performs an univariate 
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minimization. But it must be noticed that if the step is 

Pj^ , € IR (4—14) 

the sequence of tends rapidly toward 1. It therefore provides a good 
initial guess in the univariate minimization eind the line search C8in 
eventvially be abandoned in the final steps. 

The search directions are still better distributed in quasi-Newton 
methods than in first order methods. It is proven in [16] that they too 
satisfy the conjugacy property. It is also reported there, and in Dixon. 
[23]. that the method suffers an imprecise line search, and that not 
reaching the local minimum at each iteration can indeed improve the 
convergence . 

Since updating the matrix H is not a much more expensive task than 
evaluating the residual function, as it will be seen later, the method 
appears potentially better than the conjugate gradient method, whose 
advantage lies principally in the small cost associated with the 
computation of the search directions. 

4.3 Evaluation of the Computational Task. 

4.3.1 General Purpose. 

It is necessary to estimate the computational task of any 
analytical limit cycle prediction method in order to be able to comment 
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on its effectiveness. Determining the most applicable minimization 
method shall then be possible, as well as comparing alternate methods 
that can be used to study the nonlinear dynamics, such as simulation. 
Hence, the following sections will be devoted to finding rough estimates 
of the calculations undertaken in the different parts of the limit cycle 
determination process. The most general case is studied, and the 
estimates C2in be taken as upper bounds, but they shall however be 
considered as significant results. 

4.3.2 Residual Function Calculation. 

The way to compute the residual function is to first evaluate the 

residual vector as in (4-1) using (3-44), then compute the square of its 

norm eind divide it by the square of the norm of the complex amplitude 

vector. The main chore of the computation is to evaluate the forcing 

terms through the matrices Ag and Following the steps of section 

3.2.4 to build the nonlinear terms, and referring to the expression 

(3-44), the computation of one forcing term can be summarized as 

detailed in figure 4-1. Following the calculation shown there, the 

number of required operations can be estimated to the order of 

0( 21 N.n) multiplications, when all the joints are taken into account, 
J 

and not accounting for the 3 N. evaluations of describing function 

J 

terms. 

If the describing functions are evaluated on line and computed 
through a numerical integration of some available test functions, then 


pro j . on 
input axis 

3 X n op. 


comp, of sin. 
amplitude 

3 op. 



L J 

pro j . of the 
force on joint’s 
output axis 

3 X 6 op. 


Figure 4-1: Steps of the 

due to Oi 


DIDF comp, of the harmonic 

evaluation components of force 

< 3 X n op. 5 op. 



calculation of 

the force’s contribution 

on the modes’ dynamics 

18 X n op. 


Computation of the Forcing Term 


Single .Toint 
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the extra task can be evaluated on the order of 0(3 N^n) operations, 
assuming that about n points are taken to perform the numerical 
integration. The choice seems reasonable if n is about 100 as we have 
implicitly assumed. 

On the other hand, models can be available that simply describe the 
shape of the envelope of the load-displacement test curves ( [24,25]). 
In that case, closed-forms may exist that give the values of the 
describing functions taken for each joint, and the task then would be 
reduced to 0(3Nj). 

2 

The rest of the residual function evaluation requires 0(5n ) 

operations, as shown in figure 4-1 . Therefore, we can evaluate the 

2 

chore to be of the order of 0( 24 N .n + 5 n ) multiplications, which 

.1 

represents about 2.5 million operations for one function evaluation if n 
and Nj are 100 Eind 1000 respectively. Hence, it is clear that minimizing 
the number of function evaluations is of prime concern, and that the 
fastest method should be considered. 

4.3.3 Gradient Evaluation. 

The gradient can be evaluated numerically by perturbing one 
variable at a time. However two main problems arise: 

first the size of the perturbation must be chosen small enough so that 
the step taken lies in the linear domain, but also large enough so that 
rounding errors do not dominate the variations. Unfortunately, the 
choice is made blindly since the necessary information is contained in 
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the second derivative of the function. 

The second and more crucial problem is that given the number of 
operations required to evaluate the function, repeating the process 3n 
times appears unacceptable from an execution time viewpoint. 

If an analytical calculation of the gradient is undertaken, it 
would proceed as follows. 

Since: 


e( V ) e( V ) 

Y ) “ rT r rr^ 

U, U,+ U, uj 

~1 ~1 ~1 ~1 


(4-15) 


we have: 


dr 



= 2 


e( V 


de 


TTiT,a ^ 0 

U. U-+ U- U. au. 


(4-16) 


where 


de / du? 
^ 1 


Ej - Bj8Ag/au° 

-Bjayau® u[ * B^aA^/au® u| 
-Bjayau® u| - Bjayau® uj 


(4-17) 


with 


Ej = [0,..., 0, 1 0]^,the 1 component being the i^^ element. 
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-1 = 2 
1 


V ) 


de 


«r!i; 


9u 


Ir 


- 2 


■■( y ) ul'' 

u5'^u^+ u}^uj 


where : 


Ir 

de / aut = 
~ 1 


-B .a^„/au!^u^ 
j B 1 ^ 


[a-g)^i+b g.-bj 3 E. -B.raw /au^^u^ - a^ /au^^uj 

^ c 1 J p-' J P 1 '^1 q 1 

[<jB g^-bj ] E. -B.raA /au!’'uj + a^ /aul^'u, ] 

'- c 2 J q-* ~i J*" p 1 ~1 q 1 ~1 


a £ ( V )‘ de 

dr _ 2 ~ ~ 

du]^ u;*u; au. 


ilTri o. li ^ 


r( V ) u 


li 


u!f^u!f+ uj^uj 

/Vf I '^1 '^X 


where 


de / au!^= 
~ 1 


-B.aA„/au!^u_ 

J B 1 ~0 

-[(oB G„ - B.A ] E. -B.[aA /au^^U, - dJ^ /au^^u^ ] 
c 2 jq-‘~i J P i~l q i~l-‘ 


[A-cj^I+B G.-B.A ] E. -B.[aA /au!V + dJi /au!^ 

^ c 1 J p-* ~i J P 1 ~1 q 1 ~1 


li 


and finally 

o e( V de 

(-JT* r\j ^ /Vf 

_ - 2 : r 

aw u!’^u^’+ uj^uj aw 


(4-18) 


(4-19) 

(4-20) 


(4-21) 

] 


(4-22) 
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where 

0 

2 “ 5 ( 4 - 23 ) 

Assuming the residual vector has been stored, each element of the 
gradient can be computed using (4-16) through (4-23), recalling there is 
no u|^ variable. Again, the evaluation of the partial derivatives of 
the matrices Ag, and are the longest operations. Assuming the 
state of each joint is available from the computation of the function, 
figure 4-2 shows that still 0( 18 N.n) operations are required to get 
only one partial derivative of the function. The cost of computation of 
the derivatives of the describing functions can be assumed similar to 
the cost of the computation of the describing fimction, especially if a 
model is available. The analytical evaluation of the gradient might 
therefore be marginally less expensive, but accuracy is much likely to 
be gained over the numerical evaluation. 

In summary, computation of the gradient needs approximately 3n 
times the number of operations needed to compute the function, whether 
it is evaluated numerically or analytically. The latter method may be 
however faster and probably more accurate. And even though the gradient 
is expensive to get, it is still believed, based on [16] and [21], that 
gradient euid CJuasi-Newton methods are still more efficient than 


8b / dw = 



ang/aB Wx b; 


an^/aA Mx A’ 


B B + n„ B’ 


X B’ hn’ 

uH 




n A’ 
P P. 


n A’ 

“ % 


\ V°=‘l^c * \ ■" 

n A’ + n A’ 

^ Pu P % 


u is the generic variable relatively to which the derivation is 
performed. We use the following notation- = 3z/9u 

a; = pJ^(j)cos4)j. A^^= pj(j). A^= 0. B^ = 0 if u = uY 

A’ = P?(j)sin4).. A’ = 0. A* = P?(j). B’ = 0 if u = 

u J P^ q^ k'-*'-' u j 


A’ = 0. A* = 0. A’ = 0. B* = P, (j) 

u u k'*''' 


if u = u. 

J 


Fieure 4-2: Steps of the Computation of the Derivative of the 


Forcing: Term due to One Sincle Joint relative to One Variable 
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nongradient methods. 

4.3.4 Univariate Ninimlzation Sch«ne. 

The univariate search scheme is important, since multivariate 

minimization needs rather precise line searches, and since the technique 

used should provide the required precision in as few steps as possible. 

An efficient univariate minimization scheme is used by Dixon in 

[23], based on a bracketing technique and quadratic fit to improve the 

bracket. The bracketing technique consists in having the local minimum 

of a scalar function h( a } enclosed between a lower and an upper bound 

called L and M to respect Dixon’s notation. A third point a is the 

m 

current estimate of the minimum location, and the exact values of the 

function h( a ) are known at those three points. In order to improve the 

bracket, a parabola is fitted using those points and its minimum 

estimated at . If the prediction = L + E (M - L) is such that 

0.25< E <0.75, the prediction is accepted, otherwise E is replaced by 

0.25 or 0.75 in order to have the length of the bracket decrease by at 

least one qxiarter. The exact value of the scalar function h is then 
(El 

computed at a = o'- auid a new bracket is formed using the four points 

(E) 

at L, M, and * .by enclosing the one giving the smallest value 
for h between the two nearest points. 

The first bracket is obtained by taking a unit step downhill 

from the first current point a = 0. The exact value of h is computed and 
a parabola is fitted using the two points and the slope at a=0. If a 
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bracket is not formed, the process is repeated with if the 

absolute value of the estimate of the minimum is bigger than 

or with if it is smaller. 

The acceptance criterion stops the process whenever the new 
predicted value agrees with a within ela ], or the bracket has its 


m 


m ' 


size less than e a , where e is the accuracy parameter, 
m 

The univariate scalar function used to execute a line search from 
the current point 5^ in the direction is 


M “ ) = f( + “ ^) 


(4-24) 


and its derivative relative to a used to give the downhill direction and 
to fit the first parabola is 


da 




(4-25) 


where f is a generic scalar multivariate function, ^ the gradient at 


point Xj^. and Pj^ 


the line search direction. 


The accuracy of the line search given in [16] for the different 
multivariate minimization techniques is given in terms of gradient 
discrepancy: it states that a step should be accepted if 




(4-26) 
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T) can be as large as 0.9 in the BFGS method, but it has to be at least 
equal to 0.001 in the conjugate gradient method. The search termination 
test has to be converted in terms of the precision p>arameter e, since 
the computation of the gradient is too expensive. Nevertheless, the 
value of the parameter q provides a meaningful idea on what accuracy is 
required in the two methods. 

Two factors should make the line search process in the BFGS method 
faster than in the conjugate gradient method. The first factor is the 
precision of the local minimum prediction, which has to be high only in 
the Conjugate Gradient method. The second factor is the a priori 
knowledge of the location of the minimum, which is contained in the 
value of the second derivative of the function, and which is only 
derived in the BFGS method. This latter technique should therefore have 
much faster line searches. 

4.3.5 Update of the Search Direction. 

According to formula (4-7) and (4-8), the update of the search 

direction in the conjugate gradient method only requires 3x(3n) 

multiplications, which is very low. On the other hand, the update of the 

matrix and the computation of the new direction cein be estimated as 

2 

requiring at least 5x(3n) multiplications if intermediate vectors are 
being used in the update formula (4-12) and if best advantage is taken 
of the symmetry of the It is therefore obvious that if the number 

of joints is much greater than the number of retained states, the 
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computational task of updating the matrix is negligible compared to 
the task of computing the function itself. 


4.3.6 Stability Determination. 


Recalling section 3.4, the trajectory of the limit cycle is an 
ellipse, and perturbations in the ellipse plane are the only consistent 
disturbances to apply in order to check stability. Using (3-45), the 
principal axes of the ellipse have to be computed. If cj)^ is 


<b = cos ^(u^^uVii u!" 


II. II Uj II) 


(4-27) 


r i 

or in other words the single between the vectors U- and U- , and if the 
following vectors are defined with 


= 1/2 atsin ( 


2.11 u'TlI.II U^ll cosCb 

^1 ^U 

r 2 i 2 

II U,ll - II U, ll 


) + . t= 0.1 (4-28a) 


u! = U!" cos.#), 

~L ~1 1 


+ U, sinv>, 
~1 L 


1 = 0.1 


(4-28b) 


the equation of the ellipse becomes 


= uf COS0 + u5 sin0 
~0 ~1 ~2 


(4-29) 


where 0 is a dummy variable, and where U? and u5. as defined before. 


are 
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orthogonal, thus lying along the principal axes, and with norms equal 
respectively to the maximum and minimum distance from the ellipse to its 
center . 

In order to check for stability, consistent perturbations have to 
be found. Recalling (3-52), a matrix M and a vector P must be computed. 
M is the Jacobian of the residual vector e{ V ); P can be expressed as 


P = 


0 


-2w Uj - B G„ Uf 
c 2 ~1 


2w u!" - B uj 
c 2 '^1 


(4-30) 


and is closely related to the last column of the Jacobian, as can be 
seen from (4-23). 

The main task is therefore to compute the Jacobian at the limit 
cycle conditions, which takes 0( 21N .n) operations, and to invert it. 
The Jacobian is not hollow, thus its inversion should be performed using 
for example a QR decomposition algorithm. 

Once 5V = - M ^P is found, the corresp>ondir^ initial point of the 
perturbed trajectory 2|’ must be determined in the modal space. Recalling 
(3-45), which is the equation of the limit cycle trajectory, and using 
the notation as in (3-46) for the perturbed point in state space, is: 

S' = * «3!o * Re( sj' * M,) (4-31) 

This represents some projections and some other 0( n) operations. 
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The final step involves comparing the perturbed point to the ellipse, 
which is done by comparing the expression 


r ) 1 2 

II uf 11^ j L II u5 11^ J 
~1 


- 1 


(4-32) 


to 0. Hence, once the consistency condition has been found, the 
remaining operations cost only an 0( n) and they are not therefore the 
principal source of computation. 

In summary, the stability determination consists mostly in 
inverting a large 3n x 3n matrix, and this is a computationally very 
important calculation since the matrix is not hollow. 


4,4 The Siimilar Value Test. 

We are concerned in this last section about finding a more global 
type of test which could yield some information on the size of the 
different limit cycle parameters, in order to determine some region 
where they should lie, and to be able to pick ini tial guesses to start 
the limit cycle determination process. 

The system (3-37), or ainy of its alternate expressions, states that 
two complex matrices Mq and are singular, where 

Mq = I - G( 0)^g( V ) (4-33) 

Ml = I - GU^)\i Y ) (4-34) 
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The system also states some properties for the eigenvectors of Mq 
and Mj associated with the zero eigenvalue, more specifically that they 
must be respectively and . A very conservative test, however, is to 
check only for the possible singularity of the matrices without any 
regard to the eigenvectors’ directions. 

In order to perform the test, we use a general property of matrices 
of the form 

M = I - H N (4-35) 


where M, N and H are here only generic matrices. 

The property is the following ( [8]): 

The matrix M is not singular if the matrices H and N satisfy 


S 

max 


( H ) < 


1 


S ( N ) 
max ’ 


(4-35) 


where P ) represents the maximum singular value of a matrix P. 

If the test is applied to the matrices Mq and M^, which depend on the 
limit cycle parameters, it is theoretically possible to determine values 
of V and co where the test stands and where it is violated. 

If the test is true for all values of the parameters, then no limit 
cycle can occur. However, the test is very conservative and it is very 
likely that it will not stand if £iny performance has been sought in the 
design of the controller. On the other hand, the test can provide a 
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region in the variables space where a singularity can occur, and thus 
where limit cycles can occur. Hence, it provides bounds for the 
variables, and it may help choose an initial point to start the 
iteration in the limit cycle search process. 
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CHAPTER FIVE 

EXAMPLE OF ANALYTICAL DETERMINATION 
OF LIMIT CYCLES AND THEIR STABILITY 

5.1 Scope and Means of the Section. 

Examples of analytical determination of limit cycles are derived in 
this section for simple dynamical systems. They are constituted of a 
series of one or two masses connected by nonlinear springs, and they are 
under active dynamical control through full state feedback. 

The number of nonlinearities is much smaller in these examples than 
what occurs in the case of large truss structures. However, very complex 
examples are likely to complicate the understanding of the problems 
arising in the implementation of the analytical search methods. 
Furthermore, the chore associated with their derivation would be 
considerable since no actual controlled truss model is yet available, 
and extended computer capabilities would be necessary. Hence, the study 
is limited to simple problems, but this only forbids the addressing and 
the analyzing of the size issue, and still provides valuable results on 
the performances of the analytical methods. 

The simplicity of the examples limited the need for computing 
power, and an IBM PC/AT personal computer was used to implement the 
search algorithms. The software package MATLAB, [26], was extensively 
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used for its capability to handle vectors and matrices and its ability 
to manipulate real as well as complex numbers. The complex formulation 
of the limit cycle conditions, as expressed in (3-37), was therefore 
conveniently used. The software was also utilized for simulation. It 
provided graphic outputs, and the plots and tables of results relative 
to the examples derived can be found at the end of this section, 

5,2 Presentation of Tested Systems, 

5.2.1 General Architecture- 

Two systems were used, one with only one SISO nonlinear joint, the 
other one with two joints. Their general architecture is similar: they 
are constituted of masses assembled in series by the nonlinear joints 
and clamped to a fixed wall. 



Figure 5-lb: Two Joint System. 
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The masses were taken unitary. It was supposed that full state was 
available, and that each mass could be independently driven by an 
actuator. Those last hypotheses gave complete freedom in the design of 
the linear feedback controllers and thus allowed controllers that 
induced limit cycles. The joints are represented as springs in figure 
5-1, but they are nonlinear types of spring, and they were chosen to be 
identical in case 2. 


5.2.2 Idealized Joints’ Qiaracteristics. 


A simple idealized model of a joint that fits in the hypotheses 

made in Chapter 3 was chosen. Each joint is constituted of a linear 

spring of unitary stiffness attached to a simple Coulomb friction 

element with threshold F . The displacement of the Coulomb element is 

c 

limited to a bound S. The joints were supposed massless. 


11 


|-S 




•I 


F 


Figure 5-2: Model of Nonlinear .Joint. 


This type of joint can be modeled with two states, one giving the 
total deformation, x, and one giving the sliding of the friction 
element, s. The load, F, required at one end to produce the displacement 



96 


is the following 

F=K(x-s). K=1 (5-1) 

and the dynamics of the joint is 

s=x ,if|sl<S, |f 1=F and sgn( x) = sgn( F). 

(5-2) 

s = 0 

Such a model of a joint is a very simplified version of what are 
believed to be the true elements, and it more particularly displays the 
two properties we assumed in Chapter 3, which are. that it is 
"asymptotically linear", and that it has a hysteresis whose size does 
not depend on the frequency of the excitation. 

Starting from x = 0 and s = 0, the joint is truly linear and s 
stays constant as long as F does not cross the threshold F^. This 
corresponds to the linear behavior in figure 5-3a. If s is not initially 
zero, and if the load does not cross the threshold, the joint will still 
behave linearly, but around an other operating point: because of the 
Coulomb element, there is a loss of static accuracy. 

When the load crosses the threshold, the Coulomb element "gives way" and 
starts to slide. The sliding will not stop until the maximum play is 
reached, and the spring will then be streched again. When the load is 
released, the spring will start to shrink back and may be compressed, 
but the Coulomb element will not slide back toward s = 0 until the 
compression force reaches -F^: hence the load-displacement curve 

displays a hysteresis, as shown in figure 5-3b. 
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5.2.3 Describing Function for Idealized Joints. 


The use of integrators in the control design prevents, in all 
cases, any non zero steady state errors when a static disturbance is 
applied. This implies that the gain n^ need not be calculated, and that 
the gain n^ only needs to be evaluated for B = 0, since, as it will be 
shown later, the controller forces the biases to zero. In other words. 
Sinusoidal Input Describing Functions can be used here, rather than Dual 
Input Decribing Functions. 

A closed-form expression can be found for the describing functions 
of the remaining nonlinear part of the joints. If F^ is defined as: 


pNL _ pload 


- K X, 


(5-3) 


where is the total load as in (5-1), and F = K x is the linear 
spring model teiken to replace it, then we have the following expressions 
for n^: 


n.( A, 0) = n ( A, 0) + j n ( A, 0) 


Hp( A, 0) = 0 

n ( A, 0) = 0 

q 


if A < F 


(5-4a) 


«p( A. 0) = I [ -1 + f( ] 


Hq( A, 0) = 


2Fc - A) 


ifF <A<F +S (5-4b) 
c c 
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Hp( A. 0) 
Hq( A. 0) 



] 

if A > F + S 
c 


(5-4c) 


where the function f is given by: 

= -1. z < -1 


f( z) 


= [sin"^z + pi 

= +1 , z > 1 


- z^ ] . -1 < z < 1 


(5-5) 


5.2.4 Closed-loop Linear Models. 


The linear model for the joints is obviously a simple linear spring 
of unitary stiffness. F = K x. where F is specified in figure 5-2. 

The controllers implemented are multivariable PID feedback control 
systems. The integrators were found necessary for two reasons: first, 
they keep the system from settling down in a steady state that shows 
constant biases, which otherwise occured due to the loss of static 
accuracy implied by the Coulomb friction. Second, they provide the 
necessary phase lag to lead the systems to limit cycle. 

It must be noticed that the examples differ somehow from the large 
systems presented in the previous chapters, since the number of states 
is superior to the number of nonlinearities in these examples. 
Nevertheless, this only requires minor changes in the implementation of 
the method, as it will be shown later, and it does not change the way 
the problem is approached. 
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The state variables taken are the displacement x. , the rate of 
displacement and the integral of the displacement S x^dt of each 
mass. Once chosen, the feedback systems yielded respectively the 
following dynamics for the closed-loop systems: 


System 1: 


where : 


X = A, X - B K, X+ B. 

~ 1 ~ Cl 1 ~ Ji 


0 10 

0 0 1 

0-10 


B = [ 0 0 if ; B 

Jl 



(5-6) 


Kj = [ -0.5 0 -1] . 

2 ind where 

pNL ^ pload _ ^ ^ 

where is the total force to be applied at one end of the joint to 

produce a displacement x. 

The displacement of the joint d is: 


d = C, X 
1 

where : 


(5-7) 


Cj = [ 0 1 0 ] . 

The poles of the closed-loop linear system are the following: 
Xi/ 2 = -0.1761 ± j 0.8607 
X 3 = -0.6478 

The daunping ratio of the oscillatory mode is only 0.2. 
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System 2: 


where : 


X = A„ X 


B X 

Cs 2 ~ 


B. 

J2 


if 


(5-8) 



■ 0 

1 

0 

0 

0 

0 ■ 




0 

0 

1 

0 

0 

0 




0 

-2 

0 

0 

1 

0 



II 

0 

0 

0 

0 

1 

0 




0 

0 

0 

0 

0 

1 




. 0 

1 

0 

0 

-1 

0 . 



B = 
C2 

1 1 

o o 

O O 

1 0 
0 0 

0 

0 

0 

1 

1 ’^ 

[ 

0 0-1 

0 0 1 

K - 

■ 6. 

3156 

4.4587 

4.7650 

5.4444 

6.8613 

2 " 

. 5. 

3056 

6.8062 

0.5750 

6.4544 

5.5138 


0 0 0 
0 0 -1 j 


1.1250 ■ 
5.2250 _ 


The if ’s are the remaining nonlinear parts of each joint. 
The displacement of the joints, called d^’s, are: 


= X 
2 ~ 


where : 



0 1 0 0 0 0 
0-10010 


(5-9) 


The poles of the closed-loop linear system are: 
^ 1/2 - ”0-05 ± j 0.5 

^3/4 = ± J 

Xg = -3.0 

Xg = -4.0 
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The dajnping ratio of the first oscillatory mode is only 0.1. The 
damping ratio of the second mode is 0.707, and therefore, only one 
closed-loop resonant mode really appears. 

The eigenvectors associated with the first oscillatory mode are: 

Vj = [1.0000 -0.0500 -0.2475 -1.0000 0.0500 0.2475]^ (5-lOa) 

Vg = [0.0000 0.5000 -0.0500 0.0000 -0.5000 0.0500]^ (5-lOb) 


The controller was designed using e i gens t rue ture assignment 
techniques . ( [8] ) , in order to have a resonant mode where x ^ = -X 2 when 
the system is oscillating along this mode, as one can check on V, and 



The linear transfer functions can be readily derived from the state 
representations of the closed-loop systems. The transfer functions to be 
used go from each forcing term to each joint’s displacement d^ . 

In case 1, the transfer function is defined as: 



= C^(pl - Aj 


+ B KJ ^ B. 

Cl 1'' Jl 


(5-11) 


In case 2, the matrix transfer function G 2 is defined as: 




= C^CpI - 


(5-12) 
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5.3 Residual Functions. 


The set of equations (3—37) can be rebuilt for those simple 
examples by writing the following generic relations that stand for both 
cases, when the indices i, k and 1 are specified: 


Displacement and forcing terms: 

d. = B. + A.exp'^^^ 

1 1 1 


Relationship between forcing terms and displacement: 

= "b \ 

0 1 (5-14) 

d = 2 G (O)F^ + 2 G (a,)F^ . 

k ^ik k ^ik ^ 

Equating the different harmonics, the sytem of equations becomes in the 
single joint case: 


[ 1 - ngG^( 0) ] B = 0 

[ 1 - n^( A.B)Gj(ju)] A = 0 


(5-15) 
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and in the two joint case: 


[[j ?] - V ° y ] \ii 


= 0 


[ 1 0 1 _ 

[ 0 1 J 


G„( jw) X 




1 A, 


<^1 = 


(5-16) 


0 


In both applications, the bias terms B, Bj and B 2 are zero, since 
both Gj(0) = 0 and G2(0) =|^q q j . The problems therefore decouple, and 
the residual function becomes, in the single joint case 


R(A,w) = 1 1 - n^(A)G^(jw) I 
and in the two joint case 


(5-17) 


R(Aj, A 2 , 4>, w) = 


-Gp (w)n.(AJ -G (u)n.(A )+l 

^21 A A ^22 ^ ^ 


A 


1 



2 


A 


1 


2 


+ 



(5-18) 


5.4 Limit Cycle Regions. 


It is simple here to find the maximum singular value of the matrix 

where is defined as in section 3.2, since the matrix is 

diagonal, equal to diag(n.(A.)) for both examples. Therefore, S (A.) 

A. 1 max^ A' 
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is bounded by = max ( n^(A)). Figure 5-4 represents the values of 

A 

1/n^ for the values of parameters = 0.1 and S = 0.5. Figures 5-5 and 
5-6 show respectively the maximum singular value of the respective 
transfer function of case 1 eind 2. Regions where limit cycles might 
occur can therefore be determined in both examples, eind the limit cycles 
parameters will verify 
System 1: 

0.5143 rad/s < u < 1.2397 rad/s 

0.1556 < A < 1.7333 (5-19) 


System 2: 

0.3082 rad/s < w < 0.8577 rad/s 

0.0889 < Aj< 3.3778 or 0.0889 < A^< 3.3778 (5-20) 

It must be noticed, however, that the mathematics only requires one 
of the variables A^ or Ag to be within the indicated reuige, and not both 
of them simultaneously. The information is, however, still valuable, 
particularly for the frequency. 

Furthermore, the nonlinearities are only perturbation terms, and it 
can be seen that the limit cycle areas are defined by values of the 
frequency for which the input forcing terms are amplified, or in other 
words, where the linear system displays a resonant peak. 





Z<D >3Z 
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Figure 5-5: Magnitude 
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5.5 Shatpe of the Residual Function. 

In order to get some insight, the log of the residual function in 
case 1 was plotted on figure 5-7 on a 3 dimensional mesh representation. 
The scale on the frequency is logarithmic. 

\ 

Particular features of the residual function include two plateaux, 
one occur ing when A-»«>or(o-»«<>, the second one when A -» 0 or oj -» 0. 
The value of the residual on the plateaux tends toward 1. In fact, in 
both cases, either n^ or Gj(u) goes to 0, and the value of the limit can 
be directly checked from the equations. 

The importance of the linear system’s resonant mode is again 
underlined. A valley appears for w below the resonant peak frequency, 
at 0.8785 rad/s, and there is a ridge around u = TTie ridge seems 
to keep any search starting with u greater than from being 
successful, and such an attempt will end up in the plateau region. 

It appears, therefore, from the general shape of the residual 
function, that the initial conditions should be around the lower bounds 
given by the Singular Value Test, in order to reach the valley region 
during the search process aind not to be trapped in the plateau area 
where the search can drift toward points located at infinity. Caution 
should also be taken during the search process so that a step does not 
lead to frequencies well above the resonant frequency. 
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5.6 Evaluation of the Limit Cycle Prediction Method Based on The BFGS 
Algorithm. 

5.6.1 General. 

Two limit cycles were found in the first system, using the BFGS 
based limit cycle prediction method. The first one was located at: 

A = 0.2739 , <j = 0.5748 rad/s 

and the second one at: 

A = 1.1510 , u = 0.6835 rad/s. 

Simulations were performed to confirm the existence of such limit 
cycles, but only one was found, corresponding to the second set of 

conditions. The results of the simulations are shown in figure 5-8 
through 5-10 . As it can be checked there, the characteristics of the 
limit cycle found by simulation were accurately predicted by the 
analytical method: the measured frequency is at u :i£ 0.68 rad/s, and the 
measured eunplitude is A as 1.14 . 

Two limit cycles were also found in the second system, 
corresponding to the following conditions: 

Aj = 0.0698 , = 0.1391 , (}) = 3.1410 rad , w = 0.4375 rad/s, 

for the first limit cycle, and, 

Aj = 2.3441 , Ag = 4.6675 , (|) = 3.1351 rad , u = 0.4950 rad/s 


for the second one. 
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Simulations shown in figures 5-11 through 5-13 confirmed the 
existence of the second limit cycle. It can be seen in the phase plane 
representation of figures 5-11 and 5-12 that the expected elliptic form 
of the trajectory suffers some distortion that shrinks it in its middle. 
However, only the first harmonic is really present in the displacement 
of the masses, as it ceui be checked in figure 5-13. The parameters of 
the limit cycle measured from the simulation results correspond to 
Aj 2^2.30 , A 2 4.63 , 4> 3.11 . (0 0.49 rad/s 

and are therefore very close to their predicted values. 

Interestingly, the limit cycles are closely related to the resonant 
mode of the linear part the systems. The frequencies of the limit cycles 
are only slightly smaller than = 0.8785 rad/s in the first case, and 
again slightly below = 0.5025 rad/s in the second case. 

Also. calculating the displacements of the masses from the 
characteristics of the joints’ deformations in the two joint system 
case, one can find that, for the first limit cycle: 


Xj = 0.0698 cos(0.4375t) 

X 2 = -0.0693 cos(0.4375t) +0.0001 sin(0.4375t) 
and for the second limit cycle: 


Xj = 2.3441 cos(0.4950t) 

X 2 = -2.3233 cos(0.4950t) + 0.0303 sin(0.4950t) (5-22) 

Hence, it is very interesting to note that, in both cases, is 
almost equal to ’’^ 2 ' therefore, that the oscillations strongly 

resemble the oscillations of the resonant mode of the second system, 
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where always equals -x^, as one can find from the expressions (5-lOa) 
and (5-lOb). 

5.6.2 Effects of the Line Search Accuracy. 

5.6.2. 1 Tests Performed. 

Searches were undertaken, for both systems, with various values of 
the line search accuracy parameter e in order to verify that the BFGS 
method suffered rough univariate minimizations. 

With the single joint system, all searches were started at 

A=0.11 , w=0.25 rad/s, (5-23) 

and e was set successively to 10“^. 10“^, lO”^ and 10~^. The convergence 
criterion compared the value of the gradient at each iteration to 10~^^, 
and terminated the search as soon at it was below the threshold. The 
number of iterations, as well as the number of function evaluations, the 
converged point and the residual are recorded in table 5.1, and figure 
5-14 shows the residual values during the different searches. 

With the two joint system, three sets of initial conditions were 
taken. The first two sets corresponded to starting points far from a 
zero of the residual function. The first set was located at 

Aj = 0.11 : A 2 = 0.11 : (f) = IT : w = 0.3 rad/s. (5-24) 

eind the second at : 


= 0.11 : A 2 = 0.11 : (j) = IT/2; u = 0.3 rad/s. 


(5-25) 
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e was set to 0.5, 0.1 and 10 . The features of the searches are 

gathered in tables 5.2 and 5.3. The residual values during the searches 
with the first set of initial conditions are displayed in figure 5-15. 
Figure 5-16 shows the values taken by the amplitude parameters for the 
different searches, and figure 5-17 shows the frequencies. The residual 
values during the searches, started with the second set of initial 
conditions, are plotted in figure 5-18. 

The last set of initial conditions was: 

=0.6 , A 2 = 0.14 , 4)=II . cj = 0.4 rad/s (5-26) 

-1 -2 

and searches were made with e equal successively to 0.5, 10 , 10 

—5 

down to 10 . The results were gathered in table 5.4, and the residual 

values were plotted in figure 5-19. 

The convergence criterion in the two joint case was unchanged. 

5. 6. 2. 2 Discussion. 

The increase of the line search accuracy appears to have little, or 
unfavorable effect on the number of iterations. It can be seen from the 
results of the single joint model that the accuracy fixes the sequence 
of points in the search process: the first trial converges to a 

different limit cycle than the three last trials ( table 5.1). On the 
other hand, figure 5-14 shows that the sequences of points become 
similar after e has been chosen small enough. 

Results obtained on the second system show that a rough line search 
with e = 0.5 is always more successful when the search is started far 
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from a minimum (tables 5.2 and 5.3), and this same value of e only 
implies one more iteration when the search is started close to a limit 
cycle ( table 5.4). Also, figures 5-16 aind 5-17 indicate that the 
sequences of points generated by the various line search accuracies are 
different. However, figure 5-19 indicates that all the sequences of 
points converge toward a similar bound as e tends to 0. Hence, it can be 
concluded that the BFGS method might be more successful when the 
sequence of points generated in the search process is not close to the 
sequence of exact local minima which should theoretically be used: this 
conclusion meets what is reported in the literature. 

The value of e for which the method is the most efficient appears 

-3 

to depend on the problem, and is for example 10 in the first case, and 
0.5 in the second case. 

5.6.3 Influence of the Initial Conditions. 

Some assessments were already made about the values of the initial 
parameters: the singular value test provided some limits for the 
amplitudes and the frequency, and the shape of the residual function 
indicated that the search should be started with parameters equal to 
their lower bounds. 

It was also noted that the characteristics of a limit cycle were 
close to the characteristics of the first resonant mode. This can help 
fix the initial phases, as well as the relative amplitudes of the modal 


variables . 
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Information about the relative amplitudes appeared not to be really 
relevant, since searches can start far away from a minimum, and the 
values of the amplitudes can be considerably modified. Some trials, 
where the inverse of the limit cycle relative amplitude was taken 
initially, where = 0.2 and A 2 = 0.1, were even found to converge 
faster than others, started with the right relative amplitude values, 
where A ^ = 0 . 1 and A 2 = 0.2. 

On the other hand, phase initial conditions appeared to strongly 
affect the length of the search. Various trials were made, to assess the 
influence of the initial phase conditions. The searches were started 
successively with 4> = -II/2, 0 , IT/2 and IT, and the remaining variables 
were kept to 

Aj = 0.11 ; A 2 = 0.11 ; 0 ) = 0.3 rad/s. (5-27) 

The line search accuracy was set to 0.5, and the convergence threshold 
to 10 Results of the searches can be found in table 5.5, and the 
residual values were plotted in figure 5-20. 

It appears clearly from the results that the fastest search occurs 
for (j) = IT, that is, for the initial phase condition corresponding to the 
phase between the oscillations of the masses in the resonant mode. For 
4> = IT/2, the search even converges toward the second limit cycle, which, 
in term of magnitude of the defining paramaters, seems to be more 
distant from the chosen initial conditions than the first limit cycle 
( table 5.5). For the other phase initial conditions, the relative 
increase in the number of iterations compared to the one required when 
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<j) = IT, is 33%. Hence, the phase information provided by the study of the 
eigenvectors of the resonant mode appears extremely valuable, eind should 
enable the search to converge in a minimum amount of steps. 

5.7 Evaluation of the Limit Cycle Prediction Method Based on the 
Conjugate Gradient Algorithm. 

5.7.1 General. 

The evaluation of the limit cycle prediction method, using- the 
conjugate gradient algorithm, was conducted in the same manner as in the 
previous section. The exact same limit cycles were found for both 
systems. The influence of the line search accuracy, as well as the 
importance of the choice of the initial phase conditions were studied, 
as in the previous case. 

5.7.2 Effects of the Line Search Accuracy. 

5.7.2. 1 Tests Performed. 

The same evaluations as in section 5.6.2. 1. were conducted on the 
first system with the conjugate gradient based search method. Results 
were gathered in table 5.6, and the residual values during the searches 
plotted in figure 5-21. 

With the second system, only the set of initial conditions (5-27) 
was used to evaluate the performEuice of the method when started far from 
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a minimum. The line search accuracy parameter e was set to 0.5, 0.1 and 
-4 

10 . Table 5.7 and figure 5-22 show the interesting features of those 

trials. 

For the evaluation of the convergence properties of the search when 
started close to a minimum, the set of initial conditions used was: 

Aj = 0.065 : A 2 = 0.140 ;<})= IT ; u = 0.45 rad/s. 

-1 -5 

e was decreased from 10 to 10 . The results are shown in table 5.8, 

and the residual values generated during the searches plotted in figure 
5-23. 

Throughout the trials, the number of iterations was limited to 50. 
5.7. 2.2 Discussion. 

The different tests proved that the line search accuracy is a 

crucial factor in the success of the search processes based on the 

conjugate gradient method, and that, in order to be efficient, any of 

these algorithms should produce a sequence of points as close as 

possible to the sequence of exact local minima generated by the 

univariate minimizations along the line search directions. Even for this 

case with as few as four variables, the search began to be admissibly 

-4 

successful for e less than 10 , as one C2tn check in table 5.7 . Figure 

5-22 accentuates the failure of search processes using a rough line 
search. 

Another interesting result is indicated in table 5.8, and is 
confirmed by the plots of figure 5-23: there appears to be a threshold 
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—4 

for the accuracy parameter e. For e above 10 , the search fails to 

—4 

converge, or is very slow. For e below 10 . however, the sequences of 

points generated by the algorithm are close to the exact sequence of 

local minima, and the total number of iterations does not noticeably 

-4 

change. Decreasing e below 10 therefore implies only more function 
evaluations. The threshold is the saune for the two different systems 
studied, but it cannot be concluded at this point, that it would not 
change when other systems are considered. 

5.7.3 Influence of the Initial Conditions. 

The influence of the initial phase conditions was investigated as 

in section 5.6.3. The initial phase was set successively to -II/2, 0. IT/2 

and U and the remaining variables were as in (5-27). e was taken to be 
-4 

10 . The convergence test was relaxed, and a search was terminated 

whenever one of these two events occured first, the gradient became less 
than 10 , or 50 iterations were completed. Table 5.9 summarizes the 

features of the searches, and figure 5-24 shows the values of the 

residuals during the different attempts. 

(J) = IT does not lead, here, to the fastest search. Interestingly, 
when II = 0, the search converges toward the second limit cycle in the 
smallest number of iterations. In that case, the initial phase is not 
close to the value of the phase corresp>onding to the limit cycle. 
However, to ensure continuity, we allowed a point to have different 

formulations in polar coordinates. Hence, as seen in table 5.9, the 
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limit cycle found by the search started with U - 0, is defined by Aj^ = 
-2.3441, Ag = 4.6675, and (f) = -0.0065, but can otherwise be defined by 
Aj = 2.3441, Ag = 4.6675 and (j) = IT - 0.0065. Therefore, it can still be 
claimed that the fastest search occured when the initial phase condition 
was close to the phase derived from the resonant mode direction. 

Thus, with the conjugate gradient method too, the phase information 
brought by the resonant mode characteristics appears to be very 
valuable, and allows the search to converge more rapidly. 

5.8 Evaluation of the Limit Cycle Stability in the Two Tnint 

5.8.1 Stability of Limit Cycle 1. 

The state variables considered in the following are Jx^^dt, , Xj^ , 
j'x 2 dt, X 2 and Xg. where x^^ is the displacement of the first mass, and x^ 
the displacement of the second mass. 

X = [ j*x^dt , x^, Xj^ , J*X 2 dt, X 2 

The limit cycle trajectory expressed in the state space is: 

X = Uj cos(0.4375t) + uj sin(0.4375t) (5-28) 

r i 

where, and can be found from the expressions of Xj(t) and XgCt) in 
formulas (5-21), and are: 

Uj=[ 0, Aj, 0, A2sin4i/'(i), A^ + A2COs4>,-A2Sin<}) w (5-29) 

Uj=[ Aj/(J, 0,-Aj (j, (A^ + A2COs(}))/(i),-A2Sin(j),-(Aj + A2COs({)) w (5-30) 
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Using the values defining the first limit cycle, one finds: 

U[=[ 0. 6.9795 lO"^, 0, 1.7461 10”'^.-6.9318 10“^,-3.6 lO”^]^ 

Uj=[1.5953 10 0,-3.0537 10“^, -1.5844 10“^,-7.64 10~^, 3.0328 10“^]^ 

The transformation into principal axes indicated in section 4.3.6 
yields the following vectors: 

y^=[2.85 10“®,6.9795 10"^,-5.4 lO"^, 1.7178 10"f-6.9318 10~?-3.25 10“^]^ 
J^=[1.5953 10 ^-1.2 10“®-3.0537 10“?-1.5844 10”^,7.76 10“^,3.0328 10“^]^ 
The matrix M and vector P are respectively found to be 


1.0000 -1.7324 0.0000 -0.7620 
0.0000 -0.5019 -0.0698 0.5648 
0.0000 2.4534 0.0000 1.5229 
0.0000 0.9991 0.0000 -1.1292 


(5-31) 


0.5648 

0.7620 

-1.1292 

-1.5229 


(5-32) 


Perturbation parameters were found using (3-52), and expression 
(5-29) was used to generate an initial point in the state space 
corresponding to these perturbed values. The stability index p derived 
in expression (4-32) was computed, and was found to be: p = +12.04 6a. 
The positive sign of p indicates that the limit cycle is unstable . It is 
not surprising, therefore, that it could not be found through 
simulation, since infinitesimal perturbations make the system drift away 
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from the limit cycle trajectory. 


5.8.2 Stability of Limit Cycle 2. 


The p>arameters of the second limit cycle yield the following 
directions for the ellipse: 

yj=[ 0. 2.3441, 0. 6.1518 lO"^. -2.3233, -1.5074 lO"^]^ 


yj=[ 4.7355, 0, -1.1604, -4.6934, -3.0453 10^, 1.1501] 


The transformation into principal axes indicated in section 4.3.6 
yields the following vectors: 

uf=[ 3.076 10 "^, 2.3441, -7.54 10~^, 3.103 lO"^, -2.3234, -7.60 lO'^]^ 
J^=[ 4.7354, -1.523 lO"^, -1.1604, -4.6937, -1.536 10“^, 1.150l/ 

The matrix M and vector P are respectively found to be 


0.9948 -0.0013 0.0004 -6.6431 

-0.0146 -0.0072 -1.5611 46.9348 

0.0108 -0.9974 -0.0566 13.2108 

0.0287 0.0209 -7.7851 -93.9046 


(5-33) 


46.9348 

6.6431 

-93.9046 

-13.2108 


(5-34) 


As before, the stability index was computed for the new values of 
amplitudes, phase and frequency generated through (3-52). The stability 
index was found to be: p =-41.50 6a. 
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The negative sign of p indicates that the limit cycle is stable . 

The second system has. therefore, a stable, and an unstable limit 
cycle, which, recalling section 3.1.4. leads to amplitude dependent 
behaviors. These results should not. however, be surprising, since the 
system is truly linear for small amplitudes, and the linear controller 
ensures asymptotic stability in that range, and since the system is also 
"asymptotically linear". which implies that very large initial 
conditions result in bounded responses. Hence, if limit cycles exist, 
there must be one stable, since the system is stable, and there must 
also be a stability boundary that separates the truly linear behavior 
around the origin from the more nonlinear behaviors that appear for 
larger values of the state variables. 

5.9 Conclusion. 


Examples of analytical determination of limit cycles in nonlinear 
dynamical systems were derived in this chapter, and they proved that the 
analytical techniques previously derived could successfully locate the 
conditions of sustained oscillations. 

The two different techniques presented in Chapter 4 were used. The 
influence of the line search accuracy was investigated: as reported in 
[16] or [23], the BFGS method was found to perform better with a poor 
line search accuracy. The process seems, however, always successful, 
whatever level of accuracy is used in the univariate searches. On the 
other hand, the conjugate gradient method seems to require very accurate 
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line searches. It can fail to succeed if the local minima along the 
search directions are not closely approached, and thus, appears to be 
very sensitive to changes in the accuracy parameter. 

The influence of the initial conditions was also investigated. Very 
valuable results were found in the Singular Value Test. The shape of the 
residual function, for the single joint case, showed that the eimplitudes 
as well as the initial frequency should be chosen around the lower 
bounds yielded by the Singular Value Test in order for the search to be 
successful. The limit cycles were also fovmd to be very close to the 
first resonant mode of the linear parts of the systems, and this might 
provide valuable initial phase information. It was shown that starting 
with a phase in accordance with the one yielded by the resonant mode 
direction resulted in the fastest searches, with either one of the 
minimization techniques. However, the conjugate gradient method 
appeared, again, to be very sensitive to the initial phase condition, 
whereas the BFGS method was found more robust. 

Generally, the BFGS method always appeared to be noticeably more 
successful than the conjugate gradient method. When the line search 
accuracy parameters were taken to be the optimal values for both 
methods, the BFGS method consistently required less iterations. In fact, 
looking back at the results found in table 5.5, the BFGS method could 
complete a search in about 4 to 6 major iterations with the two joint 
model, whereas table 5.9 shows that, with the conjugate gradient method, 
6 major iterations correspond to the most favorable case with a relaxed 
convergence test, and that the number of major iterations is spread 
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between 6 and more than 12, depending on the initial condition chosen. 
Furthermore, more function evaluations are required in the conjugate 
gradient method, since a better line search must be performed, 
increasing accordingly the length of the process. 

Being faster, and being much more robust, the BFGS method appears 
to be. therefore, very superior to the conjugate gradient method, and 
should be retained to perform the analytical determination of limit 
cycles in dynamical systems. 
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A = 0.11 : (j = 0.25 rad/s 


Line Search 
Accuracy 
Paraumeter e 

Number 

of 

Iterations 

Number of 
function 
Evaluations 

Converged 
Point 
A ; <j 

Residual 

Function 

value 

10"^ 

12 

66 

A = 1.1510 
<0 = 0.6835 

-29 

2.6 10 

lO"^ 

9 

93 

A = 0.2739 
G) = 0.5748 

-31 

7.6 10 

10-" 

7 

84 

A = 0.2739 
G) = 0.5748 

6.7 10“^® 

10-" 

7 

95 

A = 0.2739 
G) = 0.5748 

7.4 10~^® 


Table 5.1: BFGS Method’s Performances for Different Line Search 
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Table 5.2' BFGS Method’s Performances for Different Line Search 


Accuracies in Two Mass Case with Remote Starting Point 


= 0,11 ; = 0.11 : 0 = lT/2 : g) = 0.3 rad/s 


Line Search 
Accuracy 
Parameter a 

Number 

of 

Iterations 

Number of 
function 
Evaluations 


Residual 

Function 

value 

0.5 

25 

107 

A = 2.3441 
A^= 4.6675 

4) = 3.1351 
<0 = 0.4950 

9.9 10"^ 

0.1 

26 

121 

A.= 2.3441 
A^= 4.6675 

4> = 3.1351 
<0 = 0.4950 

-29 

1.3 10 

10"^ 

25 

215 

A,= 2.3441 
A2= 4.6675 

(j) = 3.1351 
u = 0.4950 

1 

o 

00 


Table 5.3: BFCS Method’s Performances for Different Line Search 


Accuracies in Two Mass Case with Remote Starting Point 













































Line Search 
Accuracy 
Parameter e 


Number 

of 

Iterations 


Number of 
function 
Evaluations 




A-= 0.0698 
A^= 0.1391 

<j) = 3.1410 
u = 0.4375 


A,= 0.0698 
0.1391 

4> = 3.1410 
w = 0.4375 


A.= 0.0698 
A^= 0.1391 

<|> = 3.1410 
w = 0.4375 


A = 0.0698 
A^= 0.1391 

(J) = 3.1410 
u = 0.4375 


A = 0.0698 
Ag= 0.1391 

4) = 3.1410 
0) = 0.4375 


A.= 0.0698 
A2= 0.1391 

<t> = 3.1410 
w = 0.4375 


Residual 

Function 

value 


2.0 10 


2.9 10 


2.4 10 


8.1 10 


1.5 10 


1.0 10 


Table 5.4: BFGS Method's Performances for Different Line Search 


Accuracies in Two Mass Case with Accurate Start 
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0.11; A^= 0.11: (o = 0.3 rad/s: & = 0.5 


Initial 

Phase 

Condition 

Number 

of 

Iterations 

Number of 
Function 
Evaluations 

Converged 

Point 

Residual 

Function 

Value 

4> =-IT/2 

20 

89 

A^= 0.0698 
A2= 0.1391 

<|) =-3.1421 
w = 0.4375 

2.9 10“^® 

4> = 0 

20 

70 

A. =-0.0698 
A2= 0.1391 

(|) = 0.0005 
w = 0.4375 

1.6 10”^® 

= IT/2 

25 

107 

A.= 2.3441 
A2= 4.6675 

4> = 3.1351 
<i) = 0.4950 

9.9 10"^® 

<1) = IT 

15 

64 

A.= 0.0698 
A2= 0.1391 

4) = 3.1410 
w = 0.4375 

2.1 10'^® 


Table 5.5: BFGS Method Performances for Different Initial Phase 



























Line Search 
Accuracy 
Parameter e 

Number 

of 

Iterations 

Number of 
function 
Evaluations 

Converged 
Point 
A ; Cl) 

Residual 

Function 

value 

10"^ 

24 

121 

A = 1.1510 
<j = 0.6835 

-21 

2.6 10 ^ 

10"^ 

10 

82 

A = 0.2739 
w = 0.5748 

6.9 10"^5 

10-5 

9 

115 

A = 0.2739 
w = 0.5748 

1.8 10"^^ 

1 

O 

9 

141 

A = 0.2739 

1.8 lO"^"^ 


Table 5.6: Conjugate Gradient Method* s Perforngoices for Different 


Line Search Accuracies in the Single Mass Case 
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Line Search 
Accuracy 
ParEuneter e 

Number 

of 

Iterations 

Number of 
function 
Evaluations 

0.5 

50 

307 

0.1 

50 

327 

lO""^ 

42 

397 


4> = 3.1400 
0) = 0.4375 


A,= 0.0698 
0.1391 

^ = 3.1410 
<j = 0.4375 


Table 5.7 : CX? Method’s Performances for Different Line Search 


Accuracies in Two Mass Case with Remote Starting Point 

























Line Search 
Accuracy 
Parameter e 


Number 

of 

Iterations 


Number of 
function 
Evaluations 



A = 0.0698 
A^= 0.1391 

<t> = 3.1411 
w = 0.4375 


A.= 0.0698 
A^= 0.1391 

4> = 3.1410 
<0 = 0.4375 


A,= 0.0698 
A^= 0.1391 

4) = 3.1410 
w = 0.4375 


A = 0.0698 
A^= 0.1391 

4> = 3.1410 
G) = 0.4375 


A = 0.0698 
A^= 0.1391 

(]) = 3.1410 
<j = 0.4375 


Residual 

Function 

value 


1.6 10 


4.1 10 


1.8 10 


1.1 10 


1.9 10 
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A^ = 0.11; A^= 0.11: = 0.3 rad/s; e=10 


Initial 

Phase 

Condition 

Number 

of 

Iterations 

Number of 
Function 
Evaluations 

Converged 

Point 

Residual 

Function 

Value 

4) =-II/2 

50 

(grad>10“®) 

610 

A = 0.0631 
0.1341 

<|) =-2.6755 
o = 0.4432 

0.0337 

II 

O 

28 

293 

A. =-2.3441 
A2= 4.6675 

(J) =-0.0065 
(j = 0.4950 

-17 

9.9 10 

4> = IT/2 

50 

(grad>10"®) 

684 

A-= 2.3441 
A2= 4.6676 

4 > = 3.1531 
(i) = 0.4950 

8.1 10"^^ 

II 

38 

373 

A = 0.0698 
A2= 0.1391 

<t> = 3.1410 
0) = 0.4375 

-20 

3.3 10 ^ 


Table 5.9: CG Method’s Performances for Different Initial Phase 
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Xj(t) and X2(t) 



Xj(t) and X2(t) 


Figure 5-13: Two Mass case: x ^ (t), x^ft). and 
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BFGS Method: Values of the Residual during Minimization 
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ITERATION 


re 5-17: BFGS Method: Frequency during Minimization with 
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= 0-5. 0.1 and 10 . Remote Initial Point. Two Maggpg 
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'igure 5-18: BFGS Method: Values of the Residual during Minimizati 
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with 6 = 0.5, 0.1 and 10 . Second Remote Initial Point. Two Masses 
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ipjure 5 19- — BFGS Method’ Values of the Residual during Minimization 
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Fiffljre 5-21 : — OG Method: Values of the Residual during Minimization 
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Figure 5-23: CG Method: Values of the Residual during Minimization 
with fe = 10 ^j_10 10 10 10 Accurate Initial Guess 


Two Masses 
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re 5-24: CXS Method: Values of the residual during Minimization 

with (b = -IT/2. 0. 17/2. U. e = 10 Two Masses 
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CHAFIER SIX 


CXX^CLUSIOfS AND REOaOIENDATI(»S 


6.1 Conclusions on the Modeling of Laree Nonlinear Dynamirgil Systems. 

In the field of Applied Mathematics, a theoretical answer to a 
problem is considered an acceptable solution only if the required 
derivations and calculations can be actually carried out, and even 
though computers keep enlarging the capability to solve complex 
problems, the study of large dynamical systems may result in problems 
whose size makes the implementation of a solution impossible. The object 
of this work was to show a possible approach to the study of the 
dynamics of large controlled structures with nonlinear joints that would 
allow the derivation of applicable resolution techniques, especially 
concerning the determination of limit cycles. 

In Chapter 2, a general modeling framework for large structural 
systems having distributed nonlinearities was shown. The modeling 
included forming an equivalent linearized structure by replacing the 
nonlinear elements by linear ones, and performing a modal analysis on 
this altered structure, thus providing a finite state variable linear 
model. Nonlinearities were then fed back as forcii^ terms in the linear 
model. This modeling offers many advantages, among which are the use of 
existing modal analysis techniques to obtain a linear model, the global 
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way the entire system is treated, as opposed to cascading subsystems 
through nonlinear elements, the simplicity of the resulting 
representation with a linear model and a nonlinear feedback and the 
easiness to include a control feedback law. Its main disadvantage lies 
in the number of approximations that have to be made, such as taking 
approximate modes of an equivalent system, or taking only a finite 
number of them, but the accuracy can be increased at the cost of 
increasing the number of states. Hence, the modeling of the structure as 
a nonlinear feedback system appears very convenient. It can be readily 
done, requiring only a little additional work during a standard modal 
analysis, when one has to choose a linear model for each joint to 
include as the finite element equivalent stiffness. 

6.2 Analytical Determination of Limit Cycles in Large Dynamical 
Systems. 

The remainder of this work was concerned with determining 
analytically the existence of limit cycles, using the model derived in 
Qiapter 2. 

In Chapter 3, it was decided to keep only one harmonic to model the 
periodical behavior of the different state variables. This hypothesis is 
believed to be satisfactorily verified when dynamical systems such as a 
truss are considered, but should nevertheless be checked whenever the 
study of an actual structure is undertaken. 

Under the single harmonic hypothesis. Dual Input Describing 
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Functions were used to model the nonlinear load-displacement law of the 
Joints. Then, limit cycle conditions were derived, as well as stability 
conditions of the eventual limit cycles. 

In Chapter 4, an applicable resolution method was derived to detect 
whether limit cycle conditions can be met in a given system, and what 
would be the limit cycle parameters, if one is detected. The methodology 
retained was based on the minimization of a scalar residual function. 
Two minimization algorithms were pointed out in the reviewed literature: 
the BFGS method was reported to be the fastest and most successful 
minimization algorithm, and the conjugate gradient method wa& reported 
to be the most applicable for very large problems. Rough evaluations of 
the computational task made in the same chapter showed that both methods 
require a considerable number of calculations, where the computation of 
the residual function and its gradient are the most important task, 
therefore making the calculation savings of the conjugate gradient 
method irrelevant. 

The importtuice of the computation was found to come principally 
from the need to calculate the effect of each and every joint on each 
and every mode, and although the single harmonic euialysis yields a 
simple model of each joint, the order of computation was shown to be an 

0( 24N.n) to evaluate one residual in the minimization process, where N. 

V j 

is the number of distributed nonlinearities zmd n the number of modes in 
the linear model. It was also shown that the gradient computation takes 
3n times this number of operations. Furthermore, even efficient search 
techniques cannot converge in less than 3n iterations, which is the 
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number of variables that are necessary to describe a limit cycle, and 
groups of 3n. or major iterations may have to be repeated an important 
number of times. 

Examples of search for limit cycles in simple systems ( one and two 
nonlinearities) were made in Chapter 5. and they confirmed the overall 
superiority of the BEGS method. This latter approach was always found to 
be faster, and was also found very robust; that is, that the required 
number of iterations to converge only moderately depends on the initial 
conditions. 

The influence of the line search accuracy on the performance of the 
different methods was investigated. Results showed that the BEGS 
algorithm performed better with a poor determination of the minima along 
the search directions, whereas the conjugate gradient method required a 
high accuracy in order to be successful. 

The problem of finding sets of initial conditions that ensure fast 
searches was also addressed. The study showed that the first resonant 
peak of the linear part of the system plays a very important role in 
limit cycle occurrence. It led to the conclusion that the limit cycles 
must resemble the first resonant mode in terms of frequency, as well as 
phase between the oscillations, and relative amplitudes of the 
vibrations of the different parts of the structure. Most efficient sets 
of initial conditions should, accordingly, have a frequency below the 
first resonance, have small amplitudes in the neighborhood of the lower 
bounds given by the singular value test, and have initial phases derived 
from the ones found in the resonant mode directions. 
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The general conclusion of the thesis is, therefore, that the 
analytical determination of limit cycles is theoretically possible, and 
is applicable for systems of reasonable size. The computational task, 
however, rapidly grows as the number of nonlinearities N. increases. The 
number of modes on which to approximate a solution is also an important 
factor since it fixes the number of variables on which limit cycles can 
depend, and if the typical data of the CDFS I Hast experiment are taken, 
where about 1000 joints act on 100 modes, the estimates derived in 
Chapter 4 easily convince one that the size of the problem has become 
too large to be treated without further simplifications. It is.therefore 
believed that only systems smaller than the Mast experiment can be 
studied with the method presented, where every single nonlinearity is 
being considered, and alternate approaches shall be taken for larger 
problems. 

6.3 Alternate Approaches and Further Study. 

6.3.1 Complementary Studies. 

Only simple examples were derived in this study, euid they brought 
no insight on the effect of the size on the performance of the 
technique. Hence, consistently larger models should be derived 2 tnd 
studied, in order to consolidate the results emd the conclusions 
presented before about the effectiveness of analytical methods to study 
limit cycles. More precisely, it would provide data to compare the 
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performances of different minimization algorithms, and enable to 
interpolate the number of major iterations and function evaluations 
required for very large problems. It would also allow the study of 
possible improvements of the algorithms, by introducing of hybrid 
methods for example. 

6-3.2 Global and Reducc^l Joints Models. 

A field of further study could be the investigation of techniques 
that enable the concentration of the effects of groups of 
nonlinearities. This would result in feeding back a smaller number of 
global non-linearities which would result in less nonlinear forcing 
terms than in the case where each and every nonlinear element is fed 
back. The computational task would be directly affected by a decrease in 
the number of nonlinear terms fed back, and could be drastically 
reduced, 

A basic simplification that would occur in the treatement of real 
models would be to simplify the action of the forcing nonlinear terms, 
reducing it to one torque or one unique force, therefore limiting the 
action to a subspace of smaller dimension. Joints linking axial struts 
for example have only axial effects, and an elaborate routine could make 
use of this a priori information by avoiding the computation of their 
transverse components on the different modes, since those components are 
known to be zero. This could reduce considerably the computation. Only 
the action of the diagonal elements requires, in fact, a full numerical 
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treatment like the one described in figure 4.1, in order to be projected 
on the 3 dimensional space. 

Thus, it might be possible to reduce the computation in the 
analysis of limit cycles in large truss structures by simplifying the 
action of the joints, and by making best use of the geometrical 
characteristics of the structure. It might also be possible to decrease 
it by finding ways to reduce the number of nonlinear forcing terms, and 
by deriving a global representation of the effects of clusters of 
joints. The influence of those further simplifications on the solution 
accuracy should also be investigated. 

6.3.3 System Simulation. 

Simulation is almost always used in the process of analyzing 
dynamical systems. However simulation has its limitations: 

Firstly, it does not enable to forecast any trend in the behavior 
of the system, and it is theoretically necessary to try all possible 
initial conditions along with all other possible operating parameters to 
get a complete description of the system’s behavior. 

Secondly, simulation can be computationally expensive, due to the 
necessity of taking very small time steps in order to get a stable and 
accurate integration scheme. However, the cost associated with the 
analytical investigation of limit cycles may be more important than the 
one associated with the direct simulation of the system, and this latter 
approach might be worthwhile. 
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Very refined models using finite elements representation, like the 
ones shown in [2], can hardly be used to simulate very large dynamical 
systems, since they retain too many variables and cannot reasonably be 
implemented. Simple models have to be utilized, and the modeling derived 
in Chapter 2 appears to fit the purpose of simulation of large 
structures. The dynamics are approximated with the same degree of 
accuracy as when the analytical method is used, and both approaches 
should therefore yield the same solutions. The nonlinear feedback 
formulation can be readily employed in an integration scheme, and 
equation (2-22) needs very few changes to be implemented in a program. 

The main task in a simulation would be aigain to calculate the 
forcing terms and their effects on the modes. The forcing term due to 
one single joint requires the knowledge of only the joint displacement 
and rate of displacement, if a simple model as presented in [2] is 
retained. The estimate of the number of operations associated with the 
calculation can be made in the same manner as in Chapter 4, and 


indicates that 0( 8 N .n) multiplications are necessary to find the 

J 

effects of all the nonlinearities on all the modes. The rest of the 


linear effects can be evaluated with 0( 2n ) operations, which is not 


preponderant when many nonlinearities are distributed in the system. The 


number of times the evaluation is performed depends on the integration 


scheme and the number of time it is called. 


The range of integration can be expected to be long, 2 ind many 
oscillations are likely to be required in order to conclude on the 
existence of a limit cycle. A method having a small per-step truncation 


159 


error should therefore be used in order to minimize the cumulative 
error. Hamming’s method ( [28]) has such a property. The method is a 
predictor-corrector type of method. The recursive integration formula is 
stable, and it is also relatively stable, which means that the rate of 
growth of the total error during the process is always less than that of 
the solution. The use of a predictor and a corrector implies that the 
function must be evaluated twice per iteration. There is a benefit 
however in that the length of the step can be varied and in that the 
accuracy of the solution can be controlled. A fourth order Runge-Kutta 
scheme is usually chosen to start the process. The following, estimates 
of the computational chore will be based on the use of Hamming’s method. 

Data from the OOFS I mast experiment locates the lowest mode at 
0.18 Hz, corresponding to the first bending mode, and places significant 
modes up to lOOHz. Even though the controls change the natural 
frequencies, those values indicate a reasonable order of magnitude for 
the closed-loop modes. 

In order to obtain a sufficiently accurate solution, and because of 
the small per-step error of Hamming’s method, the time step could 
reasonably be estimated in the order of t = . where f is the 

maximum circular frequency. A reasonable value for the time step in the 
case of the Mast experiment should therefore be t = 0.5 milliseconds. 

The frequency of the limit cycle is believed to be around the 
lowest linear mode’s frequency, and a simulation may be run for as long 
as 50 cycles, depending on the speed with which the system tends toward 
zero, or a limit cycle. This time period corresponds to about 4.5 
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minutes in the case of the COFS I Mast, and the total number of 
evaluations of the nonlinear forcing term in that case can therefore be 
estimated to be equal to about 1.120,000. 

Comparing this figure with the computation required by the 
analytical method, it must be recalled that a single major iteration of 
the analytical process requires 3nx3nx3 times the determination of a 
nonlinear term equivalent to the nonlinear forcing term computed for the 
simulation. This means 270,000 evaluations of this equivalent nonlinear 
forcing term when the data of the Mast case are taken. Thus, the 
analytical procedure should converge in less than 5 major iterations to 
be faster than the simulation, and this is most unlikely as n keeps 
increasing. 

According to these rough estimates, the simulation approach appears 
to be an interesting alternative to the analytical determination of 
limit cycles in dynamical systems having a large number of distributed 
nonlinearities. The determining factor in both approaches is to compute 
the effect of the numerous nonlinear forcing terms on a quite important 
number of modes. Simulating the system appears to require all together 
less evaluations of the linear forcing terms. Hence, it should to be 
faster to approach the problem by means of simulation. 

The analytical method yields stable as well as unstable limit 
cycles. However, it fails to declare the absence of limit cycle, which 
can be seen more easily by simulating the system. Thus, the biggest 
advantage one method can have on the other is to be faster. Based on the 
estimates derived before, simulation was shown to be faster, but this. 
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however, might become untrue if the time step was to be reduced, and 
more importantly, if the model of the joints required more than two 
states, since it would increase the number of operations needed to 
calculate the nonlinear forcing term for every individual joint. 
Nevertheless, the effectiveness of simulation has to be considered, 
especially since alternate Integration schemes such as the ones reported 
in [ 28 ] appear to be very promising, and shall result to be very 
appropriate for simulating nonlinear feedback systems. Those techniques 
are based on a better understanding of the specificity of the problem 
they try to solve. Approaches, such as Prowler’s method, utilize as 
effectively as possible the linear character of the system. Those 
methods replace the linear continuous system by its discrete equivalent, 
thus yielding the exact homogenous response, whatever time step is used. 
The inputs to the discrete system are seunpled versions of the continuous 
inputs, and compensation is used to reduce the distortion introduced by 
the seunpling. Such integration schemes are always stable, and they are 
unusually accurate, even for large time steps. 

Therefore, limit cycles should be analyzed in significantly larger 
systems with different simulation techniques. This study would allow to 
determine the best adapted one, and it would allow to compare its 
performance with the performance of the analytical methods. Definitive 
conclusions would then be made on the real efficiency of simulation 
techniques to determine limit cycles in large actively controlled 
dynamical structures. 
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